From 44f5188ab91cb80f8a7498390b2104f6a11d8b84 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Tue, 20 Feb 2024 16:20:33 +0100 Subject: [PATCH 01/15] add libphantom-amuse.F90 --- src/utils/libphantom-amuse.F90 | 2208 ++++++++++++++++++++++++++++++++ 1 file changed, 2208 insertions(+) create mode 100644 src/utils/libphantom-amuse.F90 diff --git a/src/utils/libphantom-amuse.F90 b/src/utils/libphantom-amuse.F90 new file mode 100644 index 000000000..79fb8c89d --- /dev/null +++ b/src/utils/libphantom-amuse.F90 @@ -0,0 +1,2208 @@ +! AMUSE interface library for Phantom +! (c) 2019 - 2024 Steven Rieder + +! +! Initialize Phantom and set default parameters +! + +subroutine amuse_initialize_code_new() + use dim, only:tagline + use dim, only: maxp_hard + use memory, only:allocate_memory + use mpiutils, only:init_mpi,finalise_mpi + use initial, only:initialise,finalise,startrun,endrun + use io, only:id,master,nprocs,set_io_unit_numbers,die + use evolve, only:evol + use units, only:set_units,utime,umass,udist,unit_density + use physcon, only:gram,seconds,solarm,pc,au + implicit none + integer :: nargs + character(len=120) :: infile,logfile,evfile,dumpfile + + id = 0 + + call init_mpi(id,nprocs) + call allocate_memory(int(maxp_hard,kind=8)) + call set_io_unit_numbers() + call initialise() + call amuse_set_defaults() ! replaces reading infile + call set_units(dist=1 * au,mass=1.*solarm,G=1.) + +end subroutine amuse_initialize_code_new + +subroutine amuse_initialize_code() + use dim, only:maxp,maxp_hard,maxvxyzu + use io, only:id,nprocs,iverbose + use mpiutils, only:init_mpi + use memory, only:allocate_memory + use units, only:set_units,utime,umass,udist,unit_density + use physcon, only:gram,seconds,solarm,pc,au + use timestep, only:dtmax,dtextforce + use initial, only:initialise +#ifdef IND_TIMESTEPS + use timestep_ind, only:istepfrac,ibinnow + use part, only:ibin,ibin_old,ibin_wake +#else + use timestep, only:dtcourant,dtforce +#endif + implicit none + iverbose=5 + call init_mpi(id,nprocs) + call allocate_memory(int(maxp_hard,kind=8)) + call initialise() + call amuse_set_defaults() + call amuse_set_polyk(0.) + !print*, "maxvxyzu: ", maxvxyzu + !maxvxyzu = 4 + !call set_units(dist=50.*pc,mass=4600.*solarm,G=1.) + !call set_units(dist=1.d20,mass=1.d40,G=1.) + !call set_units(dist=1.*pc,mass=1.*solarm,G=1.) + ! call set_units(dist=1.,mass=1.,time=1.) + + !umass = 1.98892d33 * gram ! AMUSE MSun + !utime = 60 * 60 * 24 * 365.25 * 1d6 * seconds ! 10 Julian Kyr + !call set_units(time=utime,mass=umass,G=1.) + !call set_units(dist=0.1*pc,mass=1.*solarm,G=1.) + call set_units(dist=1 * au,mass=1.*solarm,G=1.) + !print*, "utime/mass/dist/udens: ", utime, umass, udist, unit_density + !dtmax = 0.01 * utime + dtextforce = 1.e-6 + +#ifdef IND_TIMESTEPS + ibin(:) = 0 + ibin_old(:) = 0 + ibin_wake(:) = 0 + istepfrac = 0 + ibinnow = 0 +#else + dtcourant = huge(dtcourant) + dtforce = huge(dtforce) +#endif + + !call amuse_initialize_wind() +end subroutine amuse_initialize_code + +subroutine amuse_set_phantom_option(name, valstring, imatch) + ! This subroutine is meant to be a replacement for read_infile + use inject, only:read_options_inject + use dust, only:read_options_dust + use growth, only:read_options_growth + use metric, only:read_options_metric + use nicil_sup, only:read_options_nicil + use dust_formation, only:read_options_dust_formation + use ptmass_radiation, only:read_options_ptmass_radiation + use eos, only:read_options_eos + use cooling, only:read_options_cooling + use ptmass, only:read_options_ptmass + use damping, only:read_options_damping + use gravwaveutils, only:read_options_gravitationalwaves + use boundary_dyn, only:read_options_boundary + implicit none + character(*), intent(inout):: name, valstring + logical:: imatch, igotall + integer:: ierr + imatch = .false. + if (.not.imatch) call read_options_inject(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_dust_formation(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_ptmass_radiation(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_dust(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_growth(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_metric(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_nicil(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_eos(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_cooling(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_damping(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_ptmass(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_gravitationalwaves(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_boundary(name,valstring,imatch,igotall,ierr) + if (.not.imatch) write(*,*) "Could not set option ", name +end subroutine amuse_set_phantom_option + +subroutine amuse_initialize_wind() + ! instead of reading a wind setup, set values here + use inject, only:read_options_inject + use dust_formation, only: read_options_dust_formation + use ptmass_radiation, only: read_options_ptmass_radiation + logical:: imatch, igotall + integer:: ierr + + call read_options_inject("sonic_type", "0", imatch, igotall, ierr) + call read_options_inject("wind_velocity", "20.", imatch, igotall, ierr) + call read_options_inject("wind_inject_radius", "2.000", imatch, igotall, ierr) ! wind injection radius (au, if 0 takes Rstar) + call read_options_inject("wind_mass_rate", "1.000E-05", imatch, igotall, ierr) ! wind mass loss rate (Msun/yr) + call read_options_inject("wind_temperature", "2500.", imatch, igotall, ierr) ! wind temperature at injection radius (K, if 0 takes Teff) + call read_options_inject("iwind_resolution", "5", imatch, igotall, ierr) ! if<>0 set number of particles on the sphere, reset particle mass + call read_options_inject("nfill_domain", "0", imatch, igotall, ierr) ! number of spheres used to set the background density profile + call read_options_inject("wind_shell_spacing", "1.000", imatch, igotall, ierr) ! desired ratio of sphere spacing to particle spacing + call read_options_inject("iboundary_spheres", "5", imatch, igotall, ierr) ! number of boundary spheres (integer) + call read_options_inject("outer_boundary", "30.", imatch, igotall, ierr) ! delete gas particles outside this radius (au) + call read_options_inject("rkill", "-1.000", imatch, igotall, ierr) ! deactivate particles outside this radius (<0 is off) + + !# options controlling dust + call read_options_dust_formation("idust_opacity", "0", imatch, igotall, ierr) ! compute dust opacity (0=off,1 (bowen), 2 (nucleation)) + + !# options controlling radiation pressure from sink particles + call read_options_ptmass_radiation("isink_radiation", "1", imatch, igotall, ierr) ! sink radiation pressure method (0=off,1=alpha,2=dust,3=alpha+dust) + call read_options_ptmass_radiation("alpha_rad", "1.000", imatch, igotall, ierr) + +end subroutine amuse_initialize_wind + +subroutine amuse_commit_parameters() +end subroutine amuse_commit_parameters + +subroutine amuse_commit_particles() + !use eos, only:ieos,init_eos + !use deriv, only:derivs + !use part, only:npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,Bevol,dBevol,& + ! dustprop,ddustprop,dustfrac,ddustevol,dens,drad,radprop,dustevol,& + ! eos_vars,metrics,pxyzu,rad + !use timestep, only:time,dtmax + !use units, only:udist,utime,umass + use initial, only:startrun + use energies, only: mtot + !use timestep, only:dtmax + implicit none + !integer :: ierr + !integer :: nactive + !double precision :: dt, dtnew_first + character(len=120) :: infile,logfile,evfile,dumpfile + !call startrun("none.in", logfile,evfile,dumpfile, .true.) + call startrun(infile,logfile,evfile,dumpfile,.true.,.true.) + print*, "total mass in code unit: ", mtot + + !dtnew_first = dtmax + !dt = 0 + !nactive = npart + print*, "COMMIT_PARTICLES, calling amuse_init_evol" + !call amuse_init_evol() + call amuse_evol(.true.) +end subroutine amuse_commit_particles + +subroutine amuse_recommit_particles() + use deriv, only:derivs + use part, only:npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,Bevol,dBevol,& + dustprop,ddustprop,dustfrac,ddustevol,dens,drad,radprop,dustevol,& + eos_vars,metrics,pxyzu,rad + use timestep, only:time,dtmax + implicit none + integer :: ierr + integer :: nactive + double precision :: dt, dtnew_first + dtnew_first = dtmax + dt = 0 + nactive = npart + call derivs(1,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,& + Bevol,dBevol,rad,drad,radprop,dustprop,ddustprop,& + dustevol,ddustevol,dustfrac,eos_vars,time,dt,dtnew_first,pxyzu,dens,metrics) + print*, "RECOMMIT_PARTICLES" +end subroutine amuse_recommit_particles + +subroutine amuse_cleanup_code() + use initial, only:endrun + implicit none + + call endrun() +end subroutine amuse_cleanup_code + +! Get npart +subroutine amuse_get_npart(npart_out, nodisabled) + use part, only:npart,xyzh + implicit none + integer, intent(out) :: npart_out + logical, intent(in) :: nodisabled + + integer :: i + + if (nodisabled) then + npart_out = 0 + do i=1,npart + if (xyzh(4,i) > 0.) then + npart_out = npart_out + 1 + endif + enddo + else + npart_out = npart + endif +end subroutine amuse_get_npart + +! Set default parameters +subroutine amuse_set_defaults() + use options, only: set_default_options,iexternalforce + implicit none + + call set_default_options() + ! A few changes to Phantom's defaults + call amuse_set_gamma(1.) + call amuse_set_ieos(1) + +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 + use evwrite, only:write_evfile,write_evlog + use energies, only:etot,totmom,angtot,mdust,np_cs_eq_0,np_e_eq_0,hdivBB_xa,& + 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,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 + 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 + + integer :: flag + character(len=256) :: infile + character(len=256) :: logfile,evfile,dumpfile + integer :: i,noutput,noutput_dtmax,nsteplast,ncount_fulldumps + real :: dtnew,dtlast,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 :: dtprint + integer :: nactive,istepfrac + integer(kind=1) :: nbinmax + logical, parameter :: dt_changed = .false. +#endif + integer :: npart_old +#ifdef INJECT_PARTICLES +#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 + + tzero = time + tlast = time + if (amuse_initialise) 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 + +#ifdef IND_TIMESTEPS + use_global_dt = .false. + istepfrac = 0 + tlast = time + write(*,*) "\n\n\n***********tlast: ", tlast, "\n\n\n" + 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 + use_global_dt = .true. + nskip = int(ntot) + nactive = npart + istepfrac = 0 ! dummy values + nbinmax = 0 + 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) + + endif !(amuse_initialise) +!end subroutine amuse_init_evol + + if (.not.amuse_initialise) then +!subroutine amuse_new_step(tlast) +! 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_dratio,check_dtmax_for_decrease +! use energies, only:etot,totmom,angtot,mdust,np_cs_eq_0,np_e_eq_0 +! use dim, only:maxvxyzu,mhd,periodic +! use fileutils, only:getnextfilename +! use options, only:nfulldump,twallmax,nmaxdumps,rhofinal1,use_dustfrac,iexternalforce,& +! icooling,ieos,ipdv_heating,ishock_heating,iresistive_heating +! 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 SORT +! use sort_particles, only:sort_part +!#endif +!#ifdef IND_TIMESTEPS +! use dim, only:maxp +! use part, only:maxphase,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 +! use timestep, only:dtdiff,C_cool +! 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 +!#ifdef MPI +! use part, only:ideadhead,shuffle_part +!#endif +!#ifdef IND_TIMESTEPS +! use part, only:twas +! use timestep_ind, only:get_dt +!#endif +! use part, only:npart,nptmass,xyzh,vxyzu,fxyzu,fext,divcurlv,massoftype, & +! xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,gravity,iboundary,npartoftype, & +! fxyz_ptmass_sinksink,ntot,poten,ndustsmall +! use quitdump, only:quit +! use ptmass, only:icreate_sinks,ptmass_create,ipart_rhomax,pt_write_sinkev +! use io_summary, only:iosum_nreal,summary_counter,summary_printout,summary_printnow +! use externalforces, only:iext_spiral +!#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 +! +! real :: dtnew,dtlast,timecheck,rhomaxold,dtmax_log_dratio +! real :: tprint,dtmaxold,dtinject +! real(kind=4) :: t1,t2,tcpu1,tcpu2,tstart,tcpustart +! real(kind=4) :: twalllast,tcpulast,twallperdump,twallused +!#ifdef IND_TIMESTEPS +! integer :: i,nalive,inbin,iamtypei +! integer(kind=1) :: nbinmaxprev +! integer(kind=8) :: nmovedtot,nalivetot +! real :: fracactive,speedup,tcheck,dtau,efficiency,tbegin +! real, intent(in) :: tlast +! real(kind=4) :: tall +! real(kind=4) :: timeperbin(0:maxbins) +! logical :: dt_changed +! integer :: iloop,npart_old +!#else +! real :: dtprint +! integer :: nactive +! logical, parameter :: dt_changed = .false. +!#endif +! logical :: use_global_dt +! integer :: j,nskip,nskipped,nevwrite_threshold,nskipped_sink,nsinkwrite_threshold +! logical :: fulldump,abortrun,abortrun_bdy,at_dump_time,writedump +! +! +! --------------------- main loop ---------------------------------------- +! + !tbegin = time + tcheck = time + npart_old = npart + + !timestepping: do while ((time < tmax).and.((nsteps < nmax) .or. (nmax < 0)).and.(rhomaxnow*rhofinal1 < 1.0)) + !write(*,*) "is istepfrac (",istepfrac,") smaller than 2**nbinmax (",2**nbinmax, ")?" + timesubstepping: do while (istepfrac < 2**nbinmax) + !write(*,*) "istepfrac (",istepfrac,") is smaller than 2**nbinmax (",2**nbinmax, "), continuing" + +#ifdef INJECT_PARTICLES + ! + ! injection of new particles into simulation + ! + !if (.not. present(flag)) then + npart_old=npart + !write(*,*) "INJECTING" + call inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,npart,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,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 + write(*,*) "dtlast: ", dtlast + + !--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 + write(*,*) "new time: ", time, tlast, istepfrac, 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) + ! We don't want to write the evfile, but we do want to calculate the energies + !call write_evfile(time,dt) + call compute_energies(time) + 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(hdivBB_xa) + 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 + !call pt_write_sinkev(nptmass,time,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,fxyz_ptmass_sinksink) +#ifdef IND_TIMESTEPS + dt_changed = .false. +#endif + write(*,*) "\n\n\n***********tlast: ", tlast, "\n\n\n" + endif +! +!--write to data file if time is right +! +! if (at_dump_time) 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 timesubstepping + +endif +!end subroutine amuse_new_step +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,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 + use timestep_ind, only:istepfrac,get_dt,nbinmax,change_nbinmax,get_newbin + use timestep, only:dt,time,dtmax + use timestep, only:C_cour,rhomaxnow + use eos, only:gamma +#endif + use timestep, only:dtextforce + use units, only:umass,udist,utime + implicit none + integer :: n, i, itype + double precision :: mass, x, y, z, vx, vy, vz, u, h + double precision :: position(3), velocity(3) +#ifdef IND_TIMESTEPS + integer(kind=1) :: nbinmaxprev + real :: dtinject + dtinject = huge(dtinject) + ! dtmax = 0.01 !TODO This is arbitrarily set. Probably bad. +#endif + ! Adding a particle of unknown temperature -> use a small cooling step + dtextforce = 1.e-6 + + itype = igas + i = npart + 1 + ! print*, "************** X position of particle ", i, x, " ***************" + position(1) = x + position(2) = y + position(3) = z + velocity(1) = vx + velocity(2) = vy + velocity(3) = vz + + if (npartoftype(itype) == 0) then + massoftype(itype) = mass + endif + call add_or_update_particle(itype,position,velocity,h, & + u,i,npart,npartoftype,xyzh,vxyzu) + !abundance(:,i) = 0. + !abundance(ih2ratio,i) = 0.5 + !abundance(iHI,i) = 1. ! assume all gas is atomic hydrogen initially + if (i == 1) then + print*, "xyz vxyz u mass = ", x, y, z, vx, vy, vz, u, mass + print*, "udist, utime, umass = ", udist, utime, umass + print*, "x vx u mass in cgs = ", x*udist, vx*udist/utime, u*udist*udist/utime/utime, mass*umass + endif + +#ifdef IND_TIMESTEPS + dtinject = C_cour * h / (gamma*(gamma-1)*u)**0.5 + nbinmaxprev = nbinmax + call get_newbin(dtinject,dtmax,nbinmax,allow_decrease=.false.) + !! not doing clever stuff: all particles go in the shortest possible bin. + !! FIXME rethink this later... + nbinmax = 3 + !if (nbinmax > nbinmaxprev) then ! update number of bins if needed + ! call change_nbinmax(nbinmax,nbinmaxprev,istepfrac,dtmax,dt) + ! print*, "nbinmax (prev), time: ", nbinmax, nbinmaxprev, time + ! print*, "npart:", npart + !endif + ! put all injected particles on shortest bin + ibin(i) = nbinmax + twas(i) = time + 0.5*get_dt(dtmax,ibin(i)) + if (i == 5) then + print*, "particle ", i, ": " + print*, x, y, z, vx, vy, vz, mass, ibin(i), twas(i) + endif +#endif + +end subroutine amuse_new_sph_particle + +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 + double precision :: mass, x, y, z, vx, vy, vz, h_smooth, u + double precision :: position(3), velocity(3) + + u = 0 + itype = idarkmatter + i = npart + 1 + position(1) = x + position(2) = y + position(3) = z + velocity(1) = vx + velocity(2) = vy + velocity(3) = vz + h_smooth = 0.1 ! TODO set this to some default + if (npartoftype(itype) == 0) then + massoftype(itype) = mass + endif + + call add_or_update_particle(itype,position,velocity,h_smooth, & + u,i,npart,npartoftype,xyzh,vxyzu) +end subroutine + +subroutine amuse_new_sink_particle(j, mass, x, y, z, vx, vy, vz, & + radius, h_smooth) + use io, only:fatal + use part, only:nptmass,maxptmass,xyzmh_ptmass,vxyz_ptmass,ihacc,ihsoft + implicit none + integer :: i, j + double precision :: mass, x, y, z, vx, vy, vz, radius, h_smooth + double precision :: position(3), velocity(3) + + nptmass = nptmass + 1 + ! Replace this with something AMUSE can handle + if (nptmass > maxptmass) call fatal('creating new sink', 'nptmass > maxptmass') + i = nptmass + j = -i + xyzmh_ptmass(:,i) = 0. + xyzmh_ptmass(1,i) = x + xyzmh_ptmass(2,i) = y + xyzmh_ptmass(3,i) = z + xyzmh_ptmass(4,i) = mass + xyzmh_ptmass(ihacc,i) = radius + xyzmh_ptmass(ihsoft,i) = h_smooth + vxyz_ptmass(1,i) = vx + vxyz_ptmass(2,i) = vy + vxyz_ptmass(3,i) = vz +end subroutine + + +subroutine amuse_delete_particle(i) + use part, only:kill_particle,xyzmh_ptmass + integer, intent(in) :: i + integer :: j + if (i == abs(i)) then + call kill_particle(i) + else + j = -i + ! Sink particles can't be killed - so we just set its mass to zero + xyzmh_ptmass(4,j) = 0 + endif +end subroutine + +subroutine amuse_get_unit_length(unit_length_out) + use units, only: udist + implicit none + double precision, intent(out) :: unit_length_out + unit_length_out = udist +end subroutine amuse_get_unit_length + +subroutine amuse_set_unit_length(unit_length_in) + use units, only: udist,utime,umass,set_units + implicit none + double precision, intent(in) :: unit_length_in + !udist = unit_length_in + !call set_units(dist=udist,time=utime,G=1.) + print*, "set_unit_length called: utime/mass/dist: ", utime, umass, udist +end subroutine amuse_set_unit_length + +subroutine amuse_get_unit_mass(unit_mass_out) + use units, only: umass + implicit none + double precision, intent(out) :: unit_mass_out + unit_mass_out = umass +end subroutine amuse_get_unit_mass + +subroutine amuse_set_unit_mass(unit_mass_in) + use units, only: umass,utime,udist,set_units + implicit none + double precision, intent(in) :: unit_mass_in + !umass = unit_mass_in + !call set_units(mass=umass,time=utime,G=1.) + print*, "set_unit_mass called: utime/mass/dist: ", utime, umass, udist +end subroutine amuse_set_unit_mass + +subroutine amuse_get_unit_time(unit_time_out) + use units, only: utime + implicit none + double precision, intent(out) :: unit_time_out + unit_time_out = utime +end subroutine amuse_get_unit_time + +subroutine amuse_set_unit_time(unit_time_in) + use units, only: utime,umass,udist,set_units + implicit none + double precision, intent(out) :: unit_time_in + !utime = unit_time_in + !call set_units(time=utime,mass=umass,G=1.) + print*, "set_unit_time called: utime/mass/dist: ", utime, umass, udist +end subroutine amuse_set_unit_time + +subroutine amuse_get_constant_solarm(solarm_out) + use physcon, only:solarm + implicit none + double precision, intent(out) :: solarm_out + solarm_out = solarm +end subroutine amuse_get_constant_solarm + +subroutine amuse_get_constant_pc(pc_out) + use physcon, only:pc + implicit none + double precision, intent(out) :: pc_out + pc_out = pc +end subroutine amuse_get_constant_pc + +subroutine amuse_get_constant_planckh(planckh_out) + use physcon, only:planckh + implicit none + double precision, intent(out) :: planckh_out + planckh_out = planckh +end subroutine amuse_get_constant_planckh + +subroutine amuse_get_potential_energy(epot_out) + use energies, only:epot + implicit none + double precision, intent(out) :: epot_out + epot_out = epot +end subroutine amuse_get_potential_energy + +subroutine amuse_get_kinetic_energy(ekin_out) + use energies, only:ekin + implicit none + double precision, intent(out) :: ekin_out + ekin_out = ekin +end subroutine amuse_get_kinetic_energy + +subroutine amuse_get_thermal_energy(etherm_out) + use energies, only:etherm + implicit none + double precision, intent(out) :: etherm_out + etherm_out = etherm +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 + dt_out = dtmax +end subroutine amuse_get_time_step + +subroutine amuse_get_number_of_sph_particles(n) + use part, only:npartoftype,igas + implicit none + integer, intent(out) :: n + logical :: nodisabled + !nodisabled = .true. + nodisabled = .false. + n = npartoftype(igas) +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 + nodisabled = .true. + call amuse_get_npart(n, nodisabled) +end subroutine amuse_get_number_of_particles + +subroutine amuse_get_time(time_out) + use timestep, only:time + implicit none + double precision, intent(out) :: time_out + time_out = time +end subroutine amuse_get_time + +subroutine amuse_get_density(i, rho) + use part, only:rhoh,iphase,massoftype,xyzh + implicit none + integer :: i + double precision :: pmassi + double precision, intent(out) :: rho + pmassi = massoftype(abs(iphase(i))) + rho = rhoh(xyzh(4,i), pmassi) +end subroutine amuse_get_density + +subroutine amuse_get_pressure(i, p) + use part, only:rhoh,iphase,massoftype,xyzh + use eos, only:ieos,equationofstate + implicit none + integer :: i, eos_type + double precision :: pmassi, ponrho, rho, spsound, x, y, z + double precision, intent(out) :: p + real :: tempi + eos_type = ieos + pmassi = massoftype(abs(iphase(i))) + call amuse_get_density(i, rho) + x = xyzh(1,i) + y = xyzh(2,i) + z = xyzh(3,i) + call equationofstate(eos_type,ponrho,spsound,rho,x,y,z,tempi) + p = ponrho * rho +end subroutine amuse_get_pressure + +subroutine amuse_get_mass(i, part_mass) + use part, only:iphase,massoftype,xyzmh_ptmass + implicit none + double precision, intent(out) :: part_mass + integer, intent(in) :: i + integer :: j + if (i == abs(i)) then + part_mass = massoftype(abs(iphase(i))) + else + j = -i + part_mass = xyzmh_ptmass(4,j) + endif +end subroutine amuse_get_mass + +subroutine amuse_get_state_gas(i, mass, x, y, z, vx, vy, vz, u, h) + implicit none + integer :: 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) + call amuse_get_velocity(i, vx, vy, vz) + call amuse_get_internal_energy(i, u) + call amuse_get_smoothing_length(i, h) +end subroutine amuse_get_state_gas + +subroutine amuse_get_state_dm(i, mass, x, y, z, vx, vy, vz) + implicit none + integer :: i + double precision :: mass, x, y, z, vx, vy, vz + call amuse_get_mass(i, mass) + call amuse_get_position(i, x, y, z) + call amuse_get_velocity(i, vx, vy, vz) + write(*,*) 'getting dm ', i +end subroutine amuse_get_state_dm + +subroutine amuse_get_state_sink(i, mass, x, y, z, vx, vy, vz, radius) + implicit none + integer :: i + double precision :: mass, x, y, z, vx, vy, vz, radius + call amuse_get_mass(i, mass) + call amuse_get_position(i, x, y, z) + call amuse_get_velocity(i, vx, vy, vz) + call amuse_get_sink_radius(i, radius) + write(*,*) 'getting sink ', i, ': radius is ', radius +end subroutine amuse_get_state_sink + +subroutine amuse_get_sink_radius(j, radius) + use part, only:xyzmh_ptmass, ihacc + implicit none + integer :: j + double precision :: radius + radius = xyzmh_ptmass(ihacc, -j) +end subroutine amuse_get_sink_radius + +subroutine amuse_get_sink_temperature(j, temperature) + use part, only:xyzmh_ptmass, iTeff + implicit none + integer, intent(in) :: j + double precision :: temperature + temperature = xyzmh_ptmass(iTeff, -j) +end subroutine + +subroutine amuse_get_sink_luminosity(j, luminosity) + use part, only:xyzmh_ptmass, iLum + implicit none + integer, intent(in) :: j + double precision :: luminosity + luminosity = xyzmh_ptmass(iLum, -j) +end subroutine + +subroutine amuse_get_position(i, x, y, z) + use part, only:xyzh,xyzmh_ptmass + implicit none + integer, intent(in) :: i + integer :: j + double precision, intent(out) :: x, y, z + if (i == abs(i)) then + x = xyzh(1, i) + y = xyzh(2, i) + z = xyzh(3, i) + else + j = -i + x = xyzmh_ptmass(1, j) + y = xyzmh_ptmass(2, j) + z = xyzmh_ptmass(3, j) + endif +end subroutine amuse_get_position + +subroutine amuse_get_velocity(i, vx, vy, vz) + use part, only:vxyzu,vxyz_ptmass + implicit none + integer, intent(in) :: i + integer :: j + double precision, intent(out) :: vx, vy, vz + if (i == abs(i)) then + vx = vxyzu(1, i) + vy = vxyzu(2, i) + vz = vxyzu(3, i) + else + j = -i + vx = vxyz_ptmass(1, j) + vy = vxyz_ptmass(2, j) + vz = vxyz_ptmass(3, j) + endif +end subroutine amuse_get_velocity + +subroutine amuse_get_acceleration(i, fx, fy, fz) + use part, only:fxyzu,fxyz_ptmass + implicit none + integer, intent(in) :: i + integer :: j + double precision, intent(out) :: fx, fy, fz + if (i == abs(i)) then + fx = fxyzu(1, i) + fy = fxyzu(2, i) + fz = fxyzu(3, i) + else + j = -i + fx = fxyz_ptmass(1, j) + fy = fxyz_ptmass(2, j) + fz = fxyz_ptmass(3, j) + endif +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 :: j + double precision, intent(out) :: h + if (i == abs(i)) then + h = xyzh(4, i) + else + j = -i + h = xyzmh_ptmass(ihsoft,j) + endif +end subroutine amuse_get_smoothing_length + +subroutine amuse_get_radius(i, radius) + implicit none + integer, intent(in) :: i + integer :: j + double precision, intent(out) :: radius + if (i == abs(i)) then + write(*,*) "not a sink" + call amuse_get_smoothing_length(i, radius) + else + write(*,*) "a sink" + j = -i + call amuse_get_sink_radius(i, radius) + endif +end subroutine amuse_get_radius + +subroutine amuse_get_internal_energy(i, u) + use dim, only:maxvxyzu + use part, only:vxyzu + implicit none + integer, intent(in) :: i + double precision, intent(out) :: u + + if (maxvxyzu >= 4) then + u = vxyzu(4, i) + else + u = 0 + endif +end subroutine + +subroutine amuse_set_hi_abundance(i, hi_abundance) + use part, only:abundance,iHI + implicit none + integer, intent(in) :: i + double precision, intent(in) :: hi_abundance + + abundance(iHI, i) = hi_abundance +end subroutine + +subroutine amuse_get_hi_abundance(i, hi_abundance) + use part, only:abundance,iHI + implicit none + integer, intent(in) :: i + double precision, intent(out) :: hi_abundance + + hi_abundance = abundance(iHI, i) +end subroutine + +subroutine amuse_set_proton_abundance(i, proton_abundance) + use part, only:abundance,iproton + implicit none + integer, intent(in) :: i + double precision, intent(in) :: proton_abundance + + abundance(iproton, i) = proton_abundance +end subroutine + +subroutine amuse_get_proton_abundance(i, proton_abundance) + use part, only:abundance,iproton + implicit none + integer, intent(in) :: i + double precision, intent(out) :: proton_abundance + + proton_abundance = abundance(iproton, i) +end subroutine + +subroutine amuse_set_electron_abundance(i, electron_abundance) + use part, only:abundance,ielectron + implicit none + integer, intent(in) :: i + double precision, intent(in) :: electron_abundance + + abundance(ielectron, i) = electron_abundance +end subroutine + +subroutine amuse_get_electron_abundance(i, electron_abundance) + use part, only:abundance,ielectron + implicit none + integer, intent(in) :: i + double precision, intent(out) :: electron_abundance + + electron_abundance = abundance(ielectron, i) +end subroutine + +subroutine amuse_set_co_abundance(i, co_abundance) + use part, only:abundance,ico + implicit none + integer, intent(in) :: i + double precision, intent(in) :: co_abundance + + abundance(ico, i) = co_abundance +end subroutine + +subroutine amuse_get_co_abundance(i, co_abundance) + use part, only:abundance,ico + implicit none + integer, intent(in) :: i + double precision, intent(out) :: co_abundance + + co_abundance = abundance(ico, i) +end subroutine + +subroutine amuse_set_h2ratio(i, h2ratio) + use part, only:abundance,iHI,ih2ratio + implicit none + integer, intent(in) :: i + double precision, intent(in) :: h2ratio + + abundance(ih2ratio, i) = h2ratio +end subroutine + +subroutine amuse_get_h2ratio(i, h2ratio) + use part, only:abundance,iHI,ih2ratio + implicit none + integer, intent(in) :: i + double precision, intent(out) :: h2ratio + + h2ratio = abundance(ih2ratio, i) +end subroutine + +subroutine amuse_set_time_step(dt_in) + use units, only:utime + use timestep, only:dtmax + implicit none + double precision, intent(in) :: dt_in + dtmax = dt_in + print*, "dtmax: ", dtmax +end subroutine + +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 :: j + if (i == abs(i)) then + massoftype(abs(iphase(i))) = part_mass + else + j = -i + xyzmh_ptmass(4,j) = part_mass + endif +end subroutine + +subroutine amuse_set_state_gas(i, mass, x, y, z, vx, vy, vz, u, h) + implicit none + integer :: 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) + call amuse_set_velocity(i, vx, vy, vz) + call amuse_set_internal_energy(i, u) + call amuse_set_smoothing_length(i, h) +end subroutine + +subroutine amuse_set_state_dm(i, mass, x, y, z, vx, vy, vz) + implicit none + integer :: i + double precision :: mass, x, y, z, vx, vy, vz + call amuse_set_mass(i, mass) + call amuse_set_position(i, x, y, z) + call amuse_set_velocity(i, vx, vy, vz) +end subroutine + +subroutine amuse_set_state_sink(i, mass, x, y, z, vx, vy, vz, radius) + implicit none + integer :: i + double precision :: mass, x, y, z, vx, vy, vz, radius + call amuse_set_mass(i, mass) + call amuse_set_position(i, x, y, z) + call amuse_set_velocity(i, vx, vy, vz) + call amuse_set_sink_radius(i, radius) +end subroutine + +subroutine amuse_set_sink_radius(j, radius) + use part, only:xyzmh_ptmass, ihacc + implicit none + integer :: j + double precision :: radius + xyzmh_ptmass(ihacc, -j) = radius +end subroutine + +subroutine amuse_set_position(i, x, y, z) + use part, only:xyzh,xyzmh_ptmass + implicit none + integer, intent(in) :: i + integer :: j + double precision, intent(in) :: x, y, z + if (i == abs(i)) then + xyzh(1, i) = x + xyzh(2, i) = y + xyzh(3, i) = z + else + j = -i + xyzmh_ptmass(1, j) = x + xyzmh_ptmass(2, j) = y + xyzmh_ptmass(3, j) = z + endif +end subroutine + +subroutine amuse_set_velocity(i, vx, vy, vz) + use part, only:vxyzu,vxyz_ptmass + implicit none + integer, intent(in) :: i + integer :: j + double precision, intent(in) :: vx, vy, vz + if (i == abs(i)) then + vxyzu(1, i) = vx + vxyzu(2, i) = vy + vxyzu(3, i) = vz + else + j = -i + vxyz_ptmass(1, j) = vx + vxyz_ptmass(2, j) = vy + vxyz_ptmass(3, j) = vz + endif +end subroutine + +subroutine amuse_set_smoothing_length(i, h) + use part, only:xyzh,xyzmh_ptmass,ihsoft + implicit none + integer, intent(in) :: i + integer :: j + double precision, intent(in) :: h + if (i == abs(i)) then + xyzh(4, i) = h + else + j = -i + xyzmh_ptmass(ihsoft, j) = h + endif +end subroutine + +subroutine amuse_set_radius(i, radius) + implicit none + integer, intent(in) :: i + integer :: j + double precision :: radius + if (i == abs(i)) then + call amuse_set_smoothing_length(i, radius) + else + j = -i + call amuse_set_sink_radius(j, radius) + endif +end subroutine + +subroutine amuse_set_sink_temperature(j, temperature) + use part, only:xyzmh_ptmass, iTeff + implicit none + integer, intent(in) :: j + double precision :: temperature + xyzmh_ptmass(iTeff, -j) = temperature +end subroutine + +subroutine amuse_set_sink_luminosity(j, luminosity) + use part, only:xyzmh_ptmass, iLum + implicit none + integer, intent(in) :: j + double precision :: luminosity + xyzmh_ptmass(iLum, -j) = luminosity +end subroutine + +subroutine amuse_set_internal_energy(i, u) + use dim, only:maxvxyzu + use part, only:vxyzu + use timestep, only:dtextforce + implicit none + integer, intent(in) :: i + double precision, intent(in) :: u + if (maxvxyzu >= 4) then + vxyzu(4, i) = u + endif + ! Changing temperature -> better use a small cooling step + dtextforce = 1.e-8 +end subroutine + +subroutine amuse_evolve_model(tmax_in) + use timestep, only:tmax, time, dt, dtmax, rhomaxnow + ! use evolvesplit, only:init_step, finalize_step +#ifdef IND_TIMESTEPS + use timestep_ind, only:istepfrac +#endif +#ifdef INJECT_PARTICLES + use inject, only:inject_particles + use part, only:npartoftype + use partinject, only:update_injected_particles +#endif + use options, only:rhofinal1 + use ptmass, only:rho_crit + use part, only:npart + use step_lf_global, only:init_step + + use part, only: xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass + implicit none + logical :: amuse_initialise + double precision, intent(in) :: tmax_in + logical :: maximum_density_reached + real :: tlast + real :: dtinject,dtlast + integer(kind=1) :: nbinmax +#ifdef INJECT_PARTICLES + integer :: npart_old +#endif +#ifndef IND_TIMESTEPS + integer :: istepfrac + istepfrac = 0 ! dummy values +#endif + dtinject = huge(dtinject) + dtlast = 0 + nbinmax = 0 + + tmax = tmax_in ! - epsilon(tmax_in) + !dtmax = (tmax - time) + + tlast = time + write(*,*) "TIMESTEPPING: evolve from ", time, " to ", tmax + timestepping: do while (time < tmax) + +#ifdef IND_TIMESTEPS + istepfrac = 0 + !print*, "*****init timestep" + !call init_step(npart, time, dtmax) +#endif + !print*, "*****new_step - Time; tmax: ", time, tmax + !call amuse_new_step(tlast) + amuse_initialise = .false. + call amuse_evol(amuse_initialise) + !print*, "*****new_step done - Time; tmax: ", time, tmax + enddo timestepping +end subroutine + +! +! Setters and getters for parameters +! + +! Setters + +subroutine amuse_set_c_courant(C_cour_in) + use timestep, only:C_cour + implicit none + double precision, intent(in) :: C_cour_in + C_cour = C_cour_in +end subroutine + +subroutine amuse_set_c_force(C_force_in) + use timestep, only:C_force + implicit none + double precision, intent(in) :: C_force_in + C_force = C_force_in +end subroutine + +subroutine amuse_set_c_cool(C_cool_in) + use timestep, only:C_cool + implicit none + double precision, intent(in) :: C_cool_in + C_cool = C_cool_in +end subroutine + +subroutine amuse_set_tolv(tolv_in) + use timestep, only:tolv + implicit none + double precision, intent(in) :: tolv_in + tolv = tolv_in +end subroutine + +subroutine amuse_set_hfact(hfact_in) + use part, only:hfact + implicit none + double precision, intent(in) :: hfact_in + hfact = hfact_in +end subroutine + +subroutine amuse_set_tolh(tolh_in) + use options, only:tolh + implicit none + double precision, intent(in) :: tolh_in + tolh = tolh_in +end subroutine + +subroutine amuse_set_tree_accuracy(tree_accuracy_in) + use kdtree, only:tree_accuracy + implicit none + double precision, intent(in) :: tree_accuracy_in + tree_accuracy = tree_accuracy_in +end subroutine + +subroutine amuse_set_alpha(alpha_in) + use options, only:alpha + implicit none + double precision, intent(in) :: alpha_in + alpha = alpha_in +end subroutine + +subroutine amuse_set_alphamax(alphamax_in) + use options, only:alphamax + implicit none + double precision, intent(in) :: alphamax_in + alphamax = alphamax_in +end subroutine + +subroutine amuse_set_beta(beta_in) + use options, only:beta + implicit none + double precision, intent(in) :: beta_in + beta = beta_in +end subroutine + +subroutine amuse_set_avdecayconst(avdecayconst_in) + use options, only:avdecayconst + implicit none + double precision, intent(in) :: avdecayconst_in + avdecayconst = avdecayconst_in +end subroutine + +subroutine amuse_set_idamp(idamp_in) + use options, only:idamp + implicit none + integer, intent(in) :: idamp_in + idamp = idamp_in +end subroutine + +subroutine amuse_set_ieos(ieos_in) + use eos, only:ieos + implicit none + integer, intent(in) :: ieos_in + ieos = ieos_in +end subroutine + +subroutine amuse_set_icooling(icooling_in) + use io, only:id,master,iprint + use options, only:icooling,iexternalforce + use part, only:h2chemistry + use chem, only:init_chem + use cooling, only:init_cooling,Tfloor + !use cooling_ism, only:init_cooling + implicit none + integer :: ierr + integer, intent(in) :: icooling_in + icooling = icooling_in + if (icooling > 0) then + Tfloor = 1 ! K + call init_cooling(id,master,iprint,ierr) + !if (h2chemistry) then + ! if (id==master) write(iprint,*) 'initialising cooling function...' + ! call init_chem() + ! call init_h2cooling() + !else + ! call init_cooling(ierr) + ! if (ierr /= 0) call fatal('initial','error initialising cooling') + !endif + endif +end subroutine + +subroutine amuse_set_polyk(polyk_in) + use eos, only:polyk + implicit none + double precision, intent(in) :: polyk_in + polyk = polyk_in +end subroutine + +subroutine amuse_set_mu(mu_in) + use eos, only:gmw + implicit none + double precision, intent(in) :: mu_in + gmw = mu_in +end subroutine + +subroutine amuse_set_rhofinal(rhofinal_in) + use options, only:rhofinal_cgs, rhofinal1 + use units, only:unit_density + implicit none + double precision, intent(in) :: rhofinal_in + rhofinal_cgs = rhofinal_in * unit_density + if (rhofinal_cgs > 0.) then + rhofinal1 = unit_density/rhofinal_cgs + else + rhofinal1 = 0.0 + endif +end subroutine + +subroutine amuse_set_rho_crit(rho_crit_in) + use units, only:unit_density + use ptmass, only:rho_crit_cgs + implicit none + double precision, intent(in) :: rho_crit_in + rho_crit_cgs = rho_crit_in * unit_density +end subroutine + +subroutine amuse_set_r_crit(r_crit_in) + use ptmass, only:r_crit + implicit none + double precision, intent(in) :: r_crit_in + r_crit = r_crit_in +end subroutine + +subroutine amuse_set_h_acc(h_acc_in) + use ptmass, only:h_acc + implicit none + double precision, intent(in) :: h_acc_in + h_acc = h_acc_in +end subroutine + +subroutine amuse_set_h_soft_sinkgas(h_soft_sinkgas_in) + use ptmass, only:h_soft_sinkgas + implicit none + double precision, intent(in) :: h_soft_sinkgas_in + h_soft_sinkgas = h_soft_sinkgas_in +end subroutine + +subroutine amuse_set_h_soft_sinksink(h_soft_sinksink_in) + use ptmass, only:h_soft_sinksink + implicit none + double precision, intent(in) :: h_soft_sinksink_in + h_soft_sinksink = h_soft_sinksink_in +end subroutine + +subroutine amuse_set_f_acc(f_acc_in) + use ptmass, only:f_acc + implicit none + double precision, intent(in) :: f_acc_in + f_acc = f_acc_in +end subroutine + +subroutine amuse_set_iexternalforce(iexternalforce_in) + use options, only:iexternalforce + implicit none + integer, intent(in) :: iexternalforce_in + iexternalforce = iexternalforce_in +end subroutine + +subroutine amuse_set_irealvisc(irealvisc_in) + use viscosity, only:irealvisc + implicit none + integer, intent(in) :: irealvisc_in + irealvisc = irealvisc_in +end subroutine + +subroutine amuse_set_shearparam(shearparam_in) + use viscosity, only:shearparam + implicit none + double precision, intent(in) :: shearparam_in + shearparam = shearparam_in +end subroutine + +subroutine amuse_set_bulkvisc(bulkvisc_in) + use viscosity, only:bulkvisc + implicit none + double precision, intent(in) :: bulkvisc_in + bulkvisc = bulkvisc_in +end subroutine + +subroutine amuse_set_gamma(gamma_in) + use eos, only:gamma + implicit none + double precision, intent(in) :: gamma_in + gamma = gamma_in +end subroutine + +subroutine amuse_set_umass(umass_in) + use units, only:umass + implicit none + double precision, intent(in) :: umass_in + umass = umass_in +end subroutine + +subroutine amuse_set_udist(udist_in) + use units, only:udist + implicit none + double precision, intent(in) :: udist_in + udist = udist_in +end subroutine + +subroutine amuse_set_utime(utime_in) + use units, only:utime + implicit none + double precision, intent(in) :: utime_in + utime = utime_in +end subroutine + +! End of Setters + +! Getters + +subroutine amuse_get_c_courant(C_cour_out) + use timestep, only:C_cour + implicit none + double precision, intent(out) :: C_cour_out + C_cour_out = C_cour +end subroutine + +subroutine amuse_get_c_force(C_force_out) + use timestep, only:C_force + implicit none + double precision, intent(out) :: C_force_out + C_force_out = C_force +end subroutine + +subroutine amuse_get_c_cool(C_cool_out) + use timestep, only:C_cool + implicit none + double precision, intent(out) :: C_cool_out + C_cool_out = C_cool +end subroutine + +subroutine amuse_get_tolv(tolv_out) + use timestep, only:tolv + implicit none + double precision, intent(out) :: tolv_out + tolv_out = tolv +end subroutine + +subroutine amuse_get_hfact(hfact_out) + use part, only:hfact + implicit none + double precision, intent(out) :: hfact_out + hfact_out = hfact +end subroutine + +subroutine amuse_get_tolh(tolh_out) + use options, only:tolh + implicit none + double precision, intent(out) :: tolh_out + tolh_out = tolh +end subroutine + +subroutine amuse_get_tree_accuracy(tree_accuracy_out) + use kdtree, only:tree_accuracy + implicit none + double precision, intent(out) :: tree_accuracy_out + tree_accuracy_out = tree_accuracy +end subroutine + +subroutine amuse_get_alpha(alpha_out) + use options, only:alpha + implicit none + double precision, intent(out) :: alpha_out + alpha_out = alpha +end subroutine + +subroutine amuse_get_alphamax(alphamax_out) + use options, only:alphamax + implicit none + double precision, intent(out) :: alphamax_out + alphamax_out = alphamax +end subroutine + +subroutine amuse_get_beta(beta_out) + use options, only:beta + implicit none + double precision, intent(out) :: beta_out + beta_out = beta +end subroutine + +subroutine amuse_get_avdecayconst(avdecayconst_out) + use options, only:avdecayconst + implicit none + double precision, intent(out) :: avdecayconst_out + avdecayconst_out = avdecayconst +end subroutine + +subroutine amuse_get_idamp(idamp_out) + use options, only:idamp + implicit none + integer, intent(out) :: idamp_out + idamp_out = idamp +end subroutine + +subroutine amuse_get_ieos(ieos_out) + use eos, only:ieos + implicit none + integer, intent(out) :: ieos_out + ieos_out = ieos +end subroutine + +subroutine amuse_get_icooling(icooling_out) + use options, only:icooling + implicit none + integer, intent(out) :: icooling_out + icooling_out = icooling +end subroutine + +subroutine amuse_get_polyk(polyk_out) + use eos, only:polyk + implicit none + double precision, intent(out) :: polyk_out + polyk_out = polyk +end subroutine + +subroutine amuse_get_mu(mu_out) + use eos, only:gmw + implicit none + double precision, intent(out) :: mu_out + mu_out = gmw +end subroutine + +subroutine amuse_get_rhofinal(rhofinal_out) + use options, only:rhofinal_cgs + use units, only:unit_density + implicit none + double precision, intent(out) :: rhofinal_out + rhofinal_out = rhofinal_cgs / unit_density +end subroutine + +subroutine amuse_get_rho_crit(rho_crit_out) + use units, only:unit_density + use ptmass, only:rho_crit_cgs + implicit none + double precision, intent(out) :: rho_crit_out + rho_crit_out = rho_crit_cgs / unit_density +end subroutine + +subroutine amuse_get_r_crit(r_crit_out) + use ptmass, only:r_crit + implicit none + double precision, intent(out) :: r_crit_out + r_crit_out = r_crit +end subroutine + +subroutine amuse_get_h_acc(h_acc_out) + use ptmass, only:h_acc + implicit none + double precision, intent(out) :: h_acc_out + h_acc_out = h_acc +end subroutine + +subroutine amuse_get_h_soft_sinkgas(h_soft_sinkgas_out) + use ptmass, only:h_soft_sinkgas + implicit none + double precision, intent(out) :: h_soft_sinkgas_out + h_soft_sinkgas_out = h_soft_sinkgas +end subroutine + +subroutine amuse_get_h_soft_sinksink(h_soft_sinksink_out) + use ptmass, only:h_soft_sinksink + implicit none + double precision, intent(out) :: h_soft_sinksink_out + h_soft_sinksink_out = h_soft_sinksink +end subroutine + +subroutine amuse_get_f_acc(f_acc_out) + use ptmass, only:f_acc + implicit none + double precision, intent(out) :: f_acc_out + f_acc_out = f_acc +end subroutine + +subroutine amuse_get_iexternalforce(iexternalforce_out) + use options, only:iexternalforce + implicit none + integer, intent(out) :: iexternalforce_out + iexternalforce_out = iexternalforce +end subroutine + +subroutine amuse_get_irealvisc(irealvisc_out) + use viscosity, only:irealvisc + implicit none + integer, intent(out) :: irealvisc_out + irealvisc_out = irealvisc +end subroutine + +subroutine amuse_get_shearparam(shearparam_out) + use viscosity, only:shearparam + implicit none + double precision, intent(out) :: shearparam_out + shearparam_out = shearparam +end subroutine + +subroutine amuse_get_bulkvisc(bulkvisc_out) + use viscosity, only:bulkvisc + implicit none + double precision, intent(out) :: bulkvisc_out + bulkvisc_out = bulkvisc +end subroutine + +subroutine amuse_get_gamma(gamma_out) + use eos, only:gamma + implicit none + double precision, intent(out) :: gamma_out + gamma_out = gamma +end subroutine + +subroutine amuse_get_umass(umass_out) + use units, only:umass + implicit none + double precision, intent(out) :: umass_out + umass_out = umass +end subroutine + +subroutine amuse_get_utime(utime_out) + use units, only:utime + implicit none + double precision, intent(out) :: utime_out + utime_out = utime +end subroutine + +subroutine amuse_get_udist(udist_out) + use units, only:udist + implicit none + double precision, intent(out) :: udist_out + udist_out = udist +end subroutine + +! End of Getters + +subroutine amuse_inject() + use timestep, only:time + use part, only:xyzh,vxyzu,npart,npartoftype + use ptmass, only:xyzmh_ptmass,vxyz_ptmass, + use inject, only:inject_particles + implicit none + real :: dtinject, dtlast + dtlast = 0 + dtinject = huge(dtinject) + + ! if npart > 0, this will also delete 'too far away particles' + ! also used to determine number of already released shells + ! time and dtlast are used together to determine 'time passed' + ! + call inject_particles(& + time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,npart,npartoftype,dtinject& + ) +end subroutine + +subroutine amuse_update_injected() + use partinject, only:update_injected_particles + use part, only:npart + use timestep_ind, only:istepfrac,nbinmax + use timestep, only:time,dtmax,dt + implicit none + integer :: npart_old + call update_injected_particles(& + npart_old,npart,istepfrac,nbinmax,time,dtmax,dt,dtinject& + ) +end subroutine From d31c87d1750801e7d45273dde19f2863c0468db9 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Tue, 20 Feb 2024 16:25:29 +0100 Subject: [PATCH 02/15] add amuse setup --- build/Makefile_setups | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/build/Makefile_setups b/build/Makefile_setups index ae562595d..80c662f38 100644 --- a/build/Makefile_setups +++ b/build/Makefile_setups @@ -1129,3 +1129,19 @@ ifeq ($(SETUP), radiativebox) RADIATION=yes PERIODIC=yes endif + +ifeq ($(SETUP), amuse) + FPPFLAGS=-DAMUSE + # H2CHEM=yes + ISOTHERMAL=no + GRAVITY=yes + IND_TIMESTEPS=yes + KNOWN_SETUP=yes + MAXPTMASS=1000 + MAXP=10000000 + # SETUPFILE= setup_wind.f90 + # ANALYSIS= utils_getneighbours.f90 utils_raytracer_all.f90 analysis_raytracer.f90 + SINK_RADIATION=yes + DUST_NUCLEATION=yes + # INJECT_PARTICLES=yes +endif From 98def797ffd683a59b4612da58d4379fa3fa7772 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Wed, 21 Feb 2024 18:49:24 +0100 Subject: [PATCH 03/15] add libphantom-amuse build targets --- build/Makefile | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/build/Makefile b/build/Makefile index ccbba721b..3a8540bf4 100644 --- a/build/Makefile +++ b/build/Makefile @@ -861,6 +861,26 @@ libphantom: libphantom.a cleanlibphantom: rm -f $(BINDIR)/libphantom.so $(LIBPHANTOM) +.PHONY: libphantom-amuse +# extra files used in libphantom that are not in OBJECTS +SRCLIB= libphantom-amuse.F90 +OBJLIB1=${SRCLIB:.f90=.o} +OBJLIB=${OBJLIB1:.F90=.o} + +LIBPHANTOMAMUSE=$(BINDIR)/libphantom-amuse.a + +libphantom-amuse.a: checksystem checkparams $(OBJECTS) $(OBJLIB) + ar rc $(LIBPHANTOMAMUSE) $(OBJLIB) $(OBJECTS) + +libphantom-amuse.so: checksystem checkparams phantom ${OBJLIB} + $(FC) -shared $(FFLAGS) $(FPPFLAGS) $(DBLFLAG) ${OBJLIB} ${OBJECTS} $(LDFLAGS) -o $(BINDIR)/libphantom-amuse.so + +libphantom-amuse: libphantom-amuse.a + +cleanlibphantom: + rm -f $(BINDIR)/libphantom-amuse.so $(LIBPHANTOMAMUSE) + + .PHONY: pyanalysis pyanalysis: libphantom.so From 9b09d37bcb80c117000767841ee4d51e5d370c67 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Wed, 21 Feb 2024 18:49:55 +0100 Subject: [PATCH 04/15] add option to prevent writing files --- src/main/options.f90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/main/options.f90 b/src/main/options.f90 index 85887742a..541b38a98 100644 --- a/src/main/options.f90 +++ b/src/main/options.f90 @@ -54,10 +54,13 @@ module options real(kind=4), public :: mcfost_keep_part character(len=80), public :: Voronoi_limits_file - ! radiation +! radiation logical, public :: exchange_radiation_energy, limit_radiation_flux, implicit_radiation logical, public :: implicit_radiation_store_drad +! library use + logical, public :: write_files + public :: set_default_options public :: ieos,idamp public :: iopacity_type From 4ba32837275c505c3d78a84d84f02603c31bf9b2 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Wed, 21 Feb 2024 18:52:05 +0100 Subject: [PATCH 05/15] don't write files when write_files is set to .false. --- src/main/evwrite.f90 | 7 +-- src/main/initial.F90 | 84 ++++++++++++++++++---------------- src/main/options.f90 | 3 ++ src/main/ptmass.F90 | 3 ++ src/main/wind.F90 | 14 +++--- src/utils/libphantom-amuse.F90 | 10 ++-- 6 files changed, 67 insertions(+), 54 deletions(-) diff --git a/src/main/evwrite.f90 b/src/main/evwrite.f90 index 8c8e5b76f..7a58de604 100644 --- a/src/main/evwrite.f90 +++ b/src/main/evwrite.f90 @@ -76,7 +76,7 @@ module evwrite subroutine init_evfile(iunit,evfile,open_file) use io, only:id,master,warning use dim, only:maxtypes,maxalpha,maxp,maxp_hard,mhd,mhd_nonideal,lightcurve - use options, only:calc_erot,ishock_heating,ipdv_heating,use_dustfrac + use options, only:calc_erot,ishock_heating,ipdv_heating,use_dustfrac,write_files use units, only:c_is_unity use part, only:igas,idust,iboundary,istar,idarkmatter,ibulge,npartoftype,ndusttypes,maxtypes use nicil, only:use_ohm,use_hall,use_ambi @@ -234,7 +234,7 @@ subroutine init_evfile(iunit,evfile,open_file) !--all threads do above, but only master writes file ! (the open_file is to prevent an .ev file from being made during the test suite) ! - if (open_file .and. id == master) then + if (write_files .and. open_file .and. id == master) then ! !--open the file for output ! @@ -351,7 +351,7 @@ subroutine write_evfile(t,dt) use energies, only:compute_energies,ev_data_update use io, only:id,master,ievfile use timestep, only:dtmax_user - use options, only:iexternalforce + use options, only:iexternalforce,write_files use externalforces,only:accretedmass1,accretedmass2 real, intent(in) :: t,dt integer :: i,j @@ -360,6 +360,7 @@ subroutine write_evfile(t,dt) call compute_energies(t) + if (.not. write_files) return if (id==master) then !--fill in additional details that are not calculated in energies.f ev_data(iev_sum,iev_dt) = dt diff --git a/src/main/initial.F90 b/src/main/initial.F90 index bbcd5a3d7..330699c03 100644 --- a/src/main/initial.F90 +++ b/src/main/initial.F90 @@ -121,7 +121,7 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread) die,fatal,id,master,nprocs,real4,warning,iverbose use externalforces, only:externalforce,initialise_externalforces,update_externalforce,& externalforce_vdependent - use options, only:iexternalforce,icooling,use_dustfrac,rhofinal1,rhofinal_cgs + use options, only:iexternalforce,icooling,use_dustfrac,rhofinal1,rhofinal_cgs,write_files use readwrite_infile, only:read_infile,write_infile use readwrite_dumps, only:read_dump,write_fulldump use part, only:npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,Bevol,dBevol,tau, tau_lucy, & @@ -551,16 +551,18 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread) #ifdef INJECT_PARTICLES call init_inject(ierr) if (ierr /= 0) call fatal('initial','error initialising particle injection') - !rename wind profile filename - inquire(file='wind_profile1D.dat',exist=iexist) - if (iexist) then - i = len(trim(dumpfile)) - if (dumpfile(i-2:i) == 'tmp') then - file1D = dumpfile(1:i-9) // '1D.dat' - else - file1D = dumpfile(1:i-5) // '1D.dat' + if (write_files) then + !rename wind profile filename + inquire(file='wind_profile1D.dat',exist=iexist) + if (iexist) then + i = len(trim(dumpfile)) + if (dumpfile(i-2:i) == 'tmp') then + file1D = dumpfile(1:i-9) // '1D.dat' + else + file1D = dumpfile(1:i-5) // '1D.dat' + endif + call rename('wind_profile1D.dat',trim(file1D)) endif - call rename('wind_profile1D.dat',trim(file1D)) endif npart_old = npart call inject_particles(time,0.,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& @@ -638,22 +640,22 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread) ! if (id==master .and. read_input_files) call write_header(2,infile,evfile,logfile,dumpfile,ntot) - call init_evfile(ievfile,evfile,.true.) + if (write_files) call init_evfile(ievfile,evfile,.true.) call write_evfile(time,dt) - if (id==master) call write_evlog(iprint) + if (write_files .and. id==master) call write_evlog(iprint) #ifdef MFLOW - call mflow_init(imflow,evfile,infile) !take evfile in input to create string.mf - call mflow_write(time, dt) + if (write_files) call mflow_init(imflow,evfile,infile) !take evfile in input to create string.mf + if (write_files) call mflow_write(time, dt) #endif #ifdef VMFLOW - call vmflow_init(ivmflow,evfile,infile) !take evfile in input to create string_v.mflowv - call vmflow_write(time, dt) + if (write_files) call vmflow_init(ivmflow,evfile,infile) !take evfile in input to create string_v.mflowv + if (write_files) call vmflow_write(time, dt) #endif #ifdef BINPOS - call binpos_init(ibinpos,evfile) !take evfile in input to create string.binpos - call binpos_write(time, dt) + if (write_files) call binpos_init(ibinpos,evfile) !take evfile in input to create string.binpos + if (write_files) call binpos_write(time, dt) #endif ! !--Determine the maximum separation of particles @@ -757,32 +759,34 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread) !--write initial conditions to output file ! if the input file ends in .tmp or .init ! - iposinit = index(dumpfile,'.init') - ipostmp = index(dumpfile,'.tmp') - if (iposinit > 0 .or. ipostmp > 0) then + if (write_files) then + iposinit = index(dumpfile,'.init') + ipostmp = index(dumpfile,'.tmp') + if (iposinit > 0 .or. ipostmp > 0) then #ifdef HDF5 - dumpfileold = trim(dumpfile)//'.h5' + dumpfileold = trim(dumpfile)//'.h5' #else - dumpfileold = dumpfile + dumpfileold = dumpfile #endif - if (iposinit > 0) then - dumpfile = trim(dumpfile(1:iposinit-1)) - else - dumpfile = trim(dumpfile(1:ipostmp-1)) - endif - call write_fulldump(time,trim(dumpfile)) - if (id==master) call write_infile(infile,logfile,evfile,trim(dumpfile),iwritein,iprint) - ! - ! delete temporary dump file - ! - call barrier_mpi() ! Ensure all procs have read temp file before deleting - inquire(file=trim(dumpfileold),exist=iexist) - if (id==master .and. iexist) then - write(iprint,"(/,a,/)") ' ---> DELETING temporary dump file '//trim(dumpfileold)//' <---' - open(unit=idisk1,file=trim(dumpfileold),status='old') - close(unit=idisk1,status='delete') + if (iposinit > 0) then + dumpfile = trim(dumpfile(1:iposinit-1)) + else + dumpfile = trim(dumpfile(1:ipostmp-1)) + endif + call write_fulldump(time,trim(dumpfile)) + if (id==master) call write_infile(infile,logfile,evfile,trim(dumpfile),iwritein,iprint) + ! + ! delete temporary dump file + ! + call barrier_mpi() ! Ensure all procs have read temp file before deleting + inquire(file=trim(dumpfileold),exist=iexist) + if (id==master .and. iexist) then + write(iprint,"(/,a,/)") ' ---> DELETING temporary dump file '//trim(dumpfileold)//' <---' + open(unit=idisk1,file=trim(dumpfileold),status='old') + close(unit=idisk1,status='delete') + endif endif - endif + endif ! (write_files) if (id==master) then call flush_warnings() diff --git a/src/main/options.f90 b/src/main/options.f90 index 541b38a98..c29da6a84 100644 --- a/src/main/options.f90 +++ b/src/main/options.f90 @@ -173,6 +173,9 @@ subroutine set_default_options ! variable composition use_var_comp = .false. + ! for use as a library + write_files = .true. + end subroutine set_default_options end module options diff --git a/src/main/ptmass.F90 b/src/main/ptmass.F90 index e36d26066..41a533edb 100644 --- a/src/main/ptmass.F90 +++ b/src/main/ptmass.F90 @@ -41,6 +41,7 @@ module ptmass use part, only:nsinkproperties,gravity,is_accretable,& ihsoft,ihacc,ispinx,ispiny,ispinz,imacc,iJ2,iReff use io, only:iscfile,iskfile,id,master + use options, only:write_files implicit none public :: init_ptmass, finish_ptmass @@ -1596,6 +1597,7 @@ subroutine init_ptmass(nptmass,logfile) integer :: i,idot character(len=150) :: filename + if (.not. write_files) return if (id /= master) return ! only do this on master thread ! !--Extract prefix & suffix @@ -1672,6 +1674,7 @@ subroutine pt_open_sinkev(num) integer :: iunit character(len=200) :: filename + if (.not. write_files) return if (id /= master) return ! only do this on master thread if (write_one_ptfile) then diff --git a/src/main/wind.F90 b/src/main/wind.F90 index 259e21e5c..5fe7a08e4 100644 --- a/src/main/wind.F90 +++ b/src/main/wind.F90 @@ -931,6 +931,7 @@ subroutine save_windprofile(r0, v0, T0, rout, tend, tcross, filename) use physcon, only:au use dust_formation, only:idust_opacity use ptmass_radiation, only:iget_tdust + use options, only:write_files real, intent(in) :: r0, v0, T0, tend, rout real, intent(out) :: tcross !time to cross the entire integration domain character(*), intent(in) :: filename @@ -955,10 +956,11 @@ subroutine save_windprofile(r0, v0, T0, rout, tend, tcross, filename) else call init_wind(r0, v0, T0, tend, state) endif - - open(unit=1337,file=filename) - call filewrite_header(1337,nwrite) - call filewrite_state(1337,nwrite, state) + if (write_files) then + open(unit=1337,file=filename) + call filewrite_header(1337,nwrite) + call filewrite_state(1337,nwrite, state) + endif eps = 0.01 iter = 0 @@ -989,7 +991,7 @@ subroutine save_windprofile(r0, v0, T0, rout, tend, tcross, filename) .or. ( abs((mu_incr -mu_base) /mu_base) > eps ) ) then writeline = writeline + 1 - call filewrite_state(1337,nwrite, state) + if (write_files) call filewrite_state(1337,nwrite, state) r_base = state%r v_base = state%v @@ -1009,7 +1011,7 @@ subroutine save_windprofile(r0, v0, T0, rout, tend, tcross, filename) else print *,'integration successful, #',iter,' iterations required, rout = ',state%r/au endif - close(1337) + if (write_files) close(1337) !stop 'save_windprofile' if (allocated(trvurho_1D)) deallocate(trvurho_1D) diff --git a/src/utils/libphantom-amuse.F90 b/src/utils/libphantom-amuse.F90 index 79fb8c89d..7a9923308 100644 --- a/src/utils/libphantom-amuse.F90 +++ b/src/utils/libphantom-amuse.F90 @@ -166,8 +166,7 @@ subroutine amuse_commit_particles() !integer :: nactive !double precision :: dt, dtnew_first character(len=120) :: infile,logfile,evfile,dumpfile - !call startrun("none.in", logfile,evfile,dumpfile, .true.) - call startrun(infile,logfile,evfile,dumpfile,.true.,.true.) + call startrun(infile,logfile,evfile,dumpfile, .true.) print*, "total mass in code unit: ", mtot !dtnew_first = dtmax @@ -227,13 +226,14 @@ end subroutine amuse_get_npart ! Set default parameters subroutine amuse_set_defaults() - use options, only: set_default_options,iexternalforce + use options, only: set_default_options,iexternalforce,write_files implicit none call set_default_options() ! A few changes to Phantom's defaults call amuse_set_gamma(1.) call amuse_set_ieos(1) + write_files=.false. end subroutine amuse_set_defaults @@ -531,7 +531,7 @@ subroutine amuse_evol(amuse_initialise) !if (.not. present(flag)) then npart_old=npart !write(*,*) "INJECTING" - call inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,npart,npartoftype,dtinject) + 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 @@ -1803,7 +1803,7 @@ subroutine amuse_set_ieos(ieos_in) subroutine amuse_set_icooling(icooling_in) use io, only:id,master,iprint use options, only:icooling,iexternalforce - use part, only:h2chemistry + use dim, only:h2chemistry use chem, only:init_chem use cooling, only:init_cooling,Tfloor !use cooling_ism, only:init_cooling From a922e6a3b35f6180dd4347f0ddc6802e8ba59c05 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Thu, 22 Feb 2024 10:44:47 +0100 Subject: [PATCH 06/15] Update Makefile_setups [skip ci] --- build/Makefile_setups | 2 -- 1 file changed, 2 deletions(-) diff --git a/build/Makefile_setups b/build/Makefile_setups index 80c662f38..e062b7d23 100644 --- a/build/Makefile_setups +++ b/build/Makefile_setups @@ -1131,7 +1131,6 @@ ifeq ($(SETUP), radiativebox) endif ifeq ($(SETUP), amuse) - FPPFLAGS=-DAMUSE # H2CHEM=yes ISOTHERMAL=no GRAVITY=yes @@ -1140,7 +1139,6 @@ ifeq ($(SETUP), amuse) MAXPTMASS=1000 MAXP=10000000 # SETUPFILE= setup_wind.f90 - # ANALYSIS= utils_getneighbours.f90 utils_raytracer_all.f90 analysis_raytracer.f90 SINK_RADIATION=yes DUST_NUCLEATION=yes # INJECT_PARTICLES=yes From 82f52232e3d62985fb8207e70d2c6221661c6379 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Mon, 15 Apr 2024 15:06:41 +0200 Subject: [PATCH 07/15] updates --- build/Makefile_setups | 7 +- src/main/timestep.f90 | 1 + src/utils/libphantom-amuse.F90 | 376 +++++++++++++++++++-------------- 3 files changed, 226 insertions(+), 158 deletions(-) diff --git a/build/Makefile_setups b/build/Makefile_setups index 80c662f38..4947ec801 100644 --- a/build/Makefile_setups +++ b/build/Makefile_setups @@ -1135,13 +1135,14 @@ ifeq ($(SETUP), amuse) # H2CHEM=yes ISOTHERMAL=no GRAVITY=yes - IND_TIMESTEPS=yes + # IND_TIMESTEPS=yes + IND_TIMESTEPS=no KNOWN_SETUP=yes MAXPTMASS=1000 MAXP=10000000 - # SETUPFILE= setup_wind.f90 + SETUPFILE= setup_wind.f90 # ANALYSIS= utils_getneighbours.f90 utils_raytracer_all.f90 analysis_raytracer.f90 SINK_RADIATION=yes DUST_NUCLEATION=yes - # INJECT_PARTICLES=yes + INJECT_PARTICLES=yes endif diff --git a/src/main/timestep.f90 b/src/main/timestep.f90 index 99bd0e172..6b09a366e 100644 --- a/src/main/timestep.f90 +++ b/src/main/timestep.f90 @@ -32,6 +32,7 @@ module timestep real :: dtmax_dratio, dtmax_max, dtmax_min, rhomaxnow real(kind=4) :: dtwallmax integer :: dtmax_ifactor,dtmax_ifactorWT + real :: dtlast public diff --git a/src/utils/libphantom-amuse.F90 b/src/utils/libphantom-amuse.F90 index 7a9923308..2e0914237 100644 --- a/src/utils/libphantom-amuse.F90 +++ b/src/utils/libphantom-amuse.F90 @@ -97,24 +97,177 @@ subroutine amuse_set_phantom_option(name, valstring, imatch) use damping, only:read_options_damping use gravwaveutils, only:read_options_gravitationalwaves use boundary_dyn, only:read_options_boundary + use timestep, only:tmax,dtmax,C_cour,C_force,C_ent,xtol,tolv,ptol,nout,nmax,& + dtwallmax,dtmax_min,dtmax_max,dtmax_dratio + use options, only:alpha,alphaB,alphamax,alphau,twallmax,tolh,rkill,rhofinal_cgs,& + psidecayfac,overcleanfac,nmaxdumps,nfulldump,limit_radiation_flux,& + ishock_heating,iresistive_heating,ireconav,ipdv_heating,iopacity_type,& + implicit_radiation,hdivbbmax_max,exchange_radiation_energy,calc_erot,& + beta,avdecayconst + use radiation_implicit, only:tol_rad,itsmax_rad,cv_type + use radiation_utils, only:kappa_cgs + use dim, only:store_dust_temperature + use viscosity, only:shearparam,irealvisc,bulkvisc + use io, only:iverbose + use part, only:hfact,ien_type implicit none character(*), intent(inout):: name, valstring logical:: imatch, igotall integer:: ierr - imatch = .false. - if (.not.imatch) call read_options_inject(name,valstring,imatch,igotall,ierr) - if (.not.imatch) call read_options_dust_formation(name,valstring,imatch,igotall,ierr) - if (.not.imatch) call read_options_ptmass_radiation(name,valstring,imatch,igotall,ierr) - if (.not.imatch) call read_options_dust(name,valstring,imatch,igotall,ierr) - if (.not.imatch) call read_options_growth(name,valstring,imatch,igotall,ierr) - if (.not.imatch) call read_options_metric(name,valstring,imatch,igotall,ierr) - if (.not.imatch) call read_options_nicil(name,valstring,imatch,igotall,ierr) - if (.not.imatch) call read_options_eos(name,valstring,imatch,igotall,ierr) - if (.not.imatch) call read_options_cooling(name,valstring,imatch,igotall,ierr) - if (.not.imatch) call read_options_damping(name,valstring,imatch,igotall,ierr) - if (.not.imatch) call read_options_ptmass(name,valstring,imatch,igotall,ierr) - if (.not.imatch) call read_options_gravitationalwaves(name,valstring,imatch,igotall,ierr) - if (.not.imatch) call read_options_boundary(name,valstring,imatch,igotall,ierr) + real :: ratio + logical :: incl_runtime2 = .false. + imatch = .true. + select case(trim(name)) + case('logfile') + ! ignored + case('dumpfile') + ! ignored + case('tmax') + read(valstring,*,iostat=ierr) tmax + case('dtmax') + read(valstring,*,iostat=ierr) dtmax + case('nmax') + read(valstring,*,iostat=ierr) nmax + case('nout') + read(valstring,*,iostat=ierr) nout + case('nmaxdumps') + read(valstring,*,iostat=ierr) nmaxdumps + case('twallmax') + read(valstring,*,iostat=ierr) twallmax + case('dtwallmax') + read(valstring,*,iostat=ierr) dtwallmax + case('iverbose') + read(valstring,*,iostat=ierr) iverbose + case('rhofinal_cgs') + read(valstring,*,iostat=ierr) rhofinal_cgs + incl_runtime2 = .true. + case('calc_erot') + read(valstring,*,iostat=ierr) calc_erot + incl_runtime2 = .true. + case('dtmax_dratio') + read(valstring,*,iostat=ierr) dtmax_dratio + incl_runtime2 = .true. + case('dtmax_max') + read(valstring,*,iostat=ierr) dtmax_max + if (dtmax_max <= 0.0) dtmax_max = dtmax + ! to prevent comparison errors from round-off + ratio = dtmax_max/dtmax + ratio = int(ratio+0.5)+0.0001 + dtmax_max = dtmax*ratio + case('dtmax_min') + read(valstring,*,iostat=ierr) dtmax_min + ! to prevent comparison errors from round-off + if (dtmax_min > epsilon(dtmax_min)) then + ratio = dtmax/dtmax_min + ratio = int(ratio+0.5)+0.0001 + dtmax_min = dtmax/ratio + endif + case('C_cour') + read(valstring,*,iostat=ierr) C_cour + case('C_force') + read(valstring,*,iostat=ierr) C_force + case('tolv') + read(valstring,*,iostat=ierr) tolv + case('C_ent') + read(valstring,*,iostat=ierr) C_ent + case('xtol') + read(valstring,*,iostat=ierr) xtol + case('ptol') + read(valstring,*,iostat=ierr) ptol + case('hfact') + read(valstring,*,iostat=ierr) hfact + case('tolh') + read(valstring,*,iostat=ierr) tolh + case('rkill') + read(valstring,*,iostat=ierr) rkill + case('nfulldump') + read(valstring,*,iostat=ierr) nfulldump + case('alpha') + read(valstring,*,iostat=ierr) alpha + case('alphamax') + read(valstring,*,iostat=ierr) alphamax + case('alphau') + read(valstring,*,iostat=ierr) alphau + case('alphaB') + read(valstring,*,iostat=ierr) alphaB + case('psidecayfac') + read(valstring,*,iostat=ierr) psidecayfac + case('overcleanfac') + read(valstring,*,iostat=ierr) overcleanfac + case('hdivbbmax_max') + read(valstring,*,iostat=ierr) hdivbbmax_max + case('beta') + read(valstring,*,iostat=ierr) beta + case('ireconav') + read(valstring,*,iostat=ierr) ireconav + case('avdecayconst') + read(valstring,*,iostat=ierr) avdecayconst + case('ipdv_heating') + read(valstring,*,iostat=ierr) ipdv_heating + case('ishock_heating') + read(valstring,*,iostat=ierr) ishock_heating + case('iresistive_heating') + read(valstring,*,iostat=ierr) iresistive_heating + case('ien_type') + read(valstring,*,iostat=ierr) ien_type + case('irealvisc') + read(valstring,*,iostat=ierr) irealvisc + case('shearparam') + read(valstring,*,iostat=ierr) shearparam + case('bulkvisc') + read(valstring,*,iostat=ierr) bulkvisc +#ifdef MCFOST + case('use_mcfost') + read(valstring,*,iostat=ierr) use_mcfost + case('Voronoi_limits_file') + read(valstring,*,iostat=ierr) Voronoi_limits_file + use_Voronoi_limits_file = .true. + case('use_mcfost_stars') + read(valstring,*,iostat=ierr) use_mcfost_stellar_parameters + case('mcfost_computes_Lacc') + read(valstring,*,iostat=ierr) mcfost_computes_Lacc + case('mcfost_uses_PdV') + read(valstring,*,iostat=ierr) mcfost_uses_PdV + case('mcfost_keep_part') + read(valstring,*,iostat=ierr) mcfost_keep_part + case('ISM') + read(valstring,*,iostat=ierr) ISM + case('mcfost_dust_subl') + read(valstring,*,iostat=ierr) mcfost_dust_subl +#endif + case('implicit_radiation') + read(valstring,*,iostat=ierr) implicit_radiation + if (implicit_radiation) store_dust_temperature = .true. + case('gas-rad_exchange') + read(valstring,*,iostat=ierr) exchange_radiation_energy + case('flux_limiter') + read(valstring,*,iostat=ierr) limit_radiation_flux + case('iopacity_type') + read(valstring,*,iostat=ierr) iopacity_type + case('cv_type') + read(valstring,*,iostat=ierr) cv_type + case('kappa_cgs') + read(valstring,*,iostat=ierr) kappa_cgs + case('tol_rad') + read(valstring,*,iostat=ierr) tol_rad + case('itsmax_rad') + read(valstring,*,iostat=ierr) itsmax_rad + case default + imatch = .false. + if (.not.imatch) call read_options_inject(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_dust_formation(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_ptmass_radiation(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_dust(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_growth(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_metric(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_nicil(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_eos(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_cooling(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_damping(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_ptmass(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_gravitationalwaves(name,valstring,imatch,igotall,ierr) + if (.not.imatch) call read_options_boundary(name,valstring,imatch,igotall,ierr) + end select if (.not.imatch) write(*,*) "Could not set option ", name end subroutine amuse_set_phantom_option @@ -243,7 +396,7 @@ subroutine amuse_evol(amuse_initialise) 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 + 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,hdivBB_xa,& compute_energies @@ -320,7 +473,7 @@ subroutine amuse_evol(amuse_initialise) character(len=256) :: infile character(len=256) :: logfile,evfile,dumpfile integer :: i,noutput,noutput_dtmax,nsteplast,ncount_fulldumps - real :: dtnew,dtlast,timecheck,rhomaxold,dtmax_log_dratio + 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 @@ -333,13 +486,14 @@ subroutine amuse_evol(amuse_initialise) 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 - integer :: npart_old #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 @@ -350,6 +504,9 @@ subroutine amuse_evol(amuse_initialise) tzero = time tlast = time + write(*,*) "AMUSE_EVOL called, amuse_initialise=", amuse_initialise + write(*,*) "dtlast", dtlast + write(*,*) "tzero, tmax: ", tzero, dtmax, tmax if (amuse_initialise) then tprint = 0. nsteps = 0 @@ -363,6 +520,7 @@ subroutine amuse_evol(amuse_initialise) call init_conservation_checks(should_conserve_energy,should_conserve_momentum,& should_conserve_angmom,should_conserve_dustmass) + noutput = 1 noutput_dtmax = 1 ncount_fulldumps = 0 @@ -378,7 +536,9 @@ subroutine amuse_evol(amuse_initialise) use_global_dt = .false. istepfrac = 0 tlast = time - write(*,*) "\n\n\n***********tlast: ", tlast, "\n\n\n" + write(*,*) + write(*,*) "***********tlast: ", tlast + write(*,*) dt = dtmax/2.**nbinmax ! use 2.0 here to allow for step too small nmovedtot = 0 tall = 0. @@ -422,107 +582,13 @@ subroutine amuse_evol(amuse_initialise) call flush(iprint) - endif !(amuse_initialise) !end subroutine amuse_init_evol - - if (.not.amuse_initialise) then -!subroutine amuse_new_step(tlast) -! 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_dratio,check_dtmax_for_decrease -! use energies, only:etot,totmom,angtot,mdust,np_cs_eq_0,np_e_eq_0 -! use dim, only:maxvxyzu,mhd,periodic -! use fileutils, only:getnextfilename -! use options, only:nfulldump,twallmax,nmaxdumps,rhofinal1,use_dustfrac,iexternalforce,& -! icooling,ieos,ipdv_heating,ishock_heating,iresistive_heating -! 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 SORT -! use sort_particles, only:sort_part -!#endif -!#ifdef IND_TIMESTEPS -! use dim, only:maxp -! use part, only:maxphase,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 -! use timestep, only:dtdiff,C_cool -! 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 -!#ifdef MPI -! use part, only:ideadhead,shuffle_part -!#endif -!#ifdef IND_TIMESTEPS -! use part, only:twas -! use timestep_ind, only:get_dt -!#endif -! use part, only:npart,nptmass,xyzh,vxyzu,fxyzu,fext,divcurlv,massoftype, & -! xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,gravity,iboundary,npartoftype, & -! fxyz_ptmass_sinksink,ntot,poten,ndustsmall -! use quitdump, only:quit -! use ptmass, only:icreate_sinks,ptmass_create,ipart_rhomax,pt_write_sinkev -! use io_summary, only:iosum_nreal,summary_counter,summary_printout,summary_printnow -! use externalforces, only:iext_spiral -!#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 -! -! real :: dtnew,dtlast,timecheck,rhomaxold,dtmax_log_dratio -! real :: tprint,dtmaxold,dtinject -! real(kind=4) :: t1,t2,tcpu1,tcpu2,tstart,tcpustart -! real(kind=4) :: twalllast,tcpulast,twallperdump,twallused -!#ifdef IND_TIMESTEPS -! integer :: i,nalive,inbin,iamtypei -! integer(kind=1) :: nbinmaxprev -! integer(kind=8) :: nmovedtot,nalivetot -! real :: fracactive,speedup,tcheck,dtau,efficiency,tbegin -! real, intent(in) :: tlast -! real(kind=4) :: tall -! real(kind=4) :: timeperbin(0:maxbins) -! logical :: dt_changed -! integer :: iloop,npart_old -!#else -! real :: dtprint -! integer :: nactive -! logical, parameter :: dt_changed = .false. -!#endif -! logical :: use_global_dt -! integer :: j,nskip,nskipped,nevwrite_threshold,nskipped_sink,nsinkwrite_threshold -! logical :: fulldump,abortrun,abortrun_bdy,at_dump_time,writedump -! -! -! --------------------- main loop ---------------------------------------- -! - !tbegin = time - tcheck = time - npart_old = npart - - !timestepping: do while ((time < tmax).and.((nsteps < nmax) .or. (nmax < 0)).and.(rhomaxnow*rhofinal1 < 1.0)) - !write(*,*) "is istepfrac (",istepfrac,") smaller than 2**nbinmax (",2**nbinmax, ")?" - timesubstepping: do while (istepfrac < 2**nbinmax) - !write(*,*) "istepfrac (",istepfrac,") is smaller than 2**nbinmax (",2**nbinmax, "), continuing" + else ! i.e.: .not. amuse_initialise + !tcheck = time + !npart_old = npart + timesubstepping: do while ((time + 0.01 * dtmax < tmax).and.((nsteps < nmax) .or. (nmax < 0)).and.(rhomaxnow*rhofinal1 < 1.0)) + !dtmax = min(dtmax, tmax-time) ! tried this, doesn't work. + !timesubstepping: do while (istepfrac < 2**nbinmax) #ifdef INJECT_PARTICLES ! @@ -608,7 +674,6 @@ subroutine amuse_evol(amuse_initialise) endif dtlast = dt - write(*,*) "dtlast: ", dtlast !--timings for step call call get_timings(t2,tcpu2) @@ -620,7 +685,6 @@ subroutine amuse_evol(amuse_initialise) !--update time in way that is free of round-off errors time = tlast + istepfrac/real(2**nbinmaxprev)*dtmaxold - write(*,*) "new time: ", time, tlast, istepfrac, nbinmaxprev, dtmaxold !--print efficiency of partial timestep if (id==master) call print_dtind_efficiency(iverbose,nalivetot,nactivetot,tall,t2-t1,1) @@ -733,7 +797,9 @@ subroutine amuse_evol(amuse_initialise) #ifdef IND_TIMESTEPS dt_changed = .false. #endif - write(*,*) "\n\n\n***********tlast: ", tlast, "\n\n\n" + write(*,*) + write(*,*) "***********tlast, dtlast: ", tlast, dtlast + write(*,*) endif ! !--write to data file if time is right @@ -1649,7 +1715,7 @@ subroutine amuse_set_internal_energy(i, u) end subroutine subroutine amuse_evolve_model(tmax_in) - use timestep, only:tmax, time, dt, dtmax, rhomaxnow + use timestep, only:tmax, time, dt, dtmax, rhomaxnow, dtlast ! use evolvesplit, only:init_step, finalize_step #ifdef IND_TIMESTEPS use timestep_ind, only:istepfrac @@ -1670,7 +1736,7 @@ subroutine amuse_evolve_model(tmax_in) double precision, intent(in) :: tmax_in logical :: maximum_density_reached real :: tlast - real :: dtinject,dtlast + real :: dtinject integer(kind=1) :: nbinmax #ifdef INJECT_PARTICLES integer :: npart_old @@ -1680,7 +1746,7 @@ subroutine amuse_evolve_model(tmax_in) istepfrac = 0 ! dummy values #endif dtinject = huge(dtinject) - dtlast = 0 + ! dtlast = 0 nbinmax = 0 tmax = tmax_in ! - epsilon(tmax_in) @@ -1688,7 +1754,7 @@ subroutine amuse_evolve_model(tmax_in) tlast = time write(*,*) "TIMESTEPPING: evolve from ", time, " to ", tmax - timestepping: do while (time < tmax) + timestepping: do while (time+0.01 * dtmax < tmax) #ifdef IND_TIMESTEPS istepfrac = 0 @@ -2176,33 +2242,33 @@ subroutine amuse_get_udist(udist_out) ! End of Getters -subroutine amuse_inject() - use timestep, only:time - use part, only:xyzh,vxyzu,npart,npartoftype - use ptmass, only:xyzmh_ptmass,vxyz_ptmass, - use inject, only:inject_particles - implicit none - real :: dtinject, dtlast - dtlast = 0 - dtinject = huge(dtinject) - - ! if npart > 0, this will also delete 'too far away particles' - ! also used to determine number of already released shells - ! time and dtlast are used together to determine 'time passed' - ! - call inject_particles(& - time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,npart,npartoftype,dtinject& - ) -end subroutine - -subroutine amuse_update_injected() - use partinject, only:update_injected_particles - use part, only:npart - use timestep_ind, only:istepfrac,nbinmax - use timestep, only:time,dtmax,dt - implicit none - integer :: npart_old - call update_injected_particles(& - npart_old,npart,istepfrac,nbinmax,time,dtmax,dt,dtinject& - ) -end subroutine +!subroutine amuse_inject() +! use timestep, only:time +! use part, only:xyzh,vxyzu,npart,npartoftype +! use ptmass, only:xyzmh_ptmass,vxyz_ptmass, +! use inject, only:inject_particles +! implicit none +! real :: dtinject, dtlast +! dtlast = 0 +! dtinject = huge(dtinject) +! +! ! if npart > 0, this will also delete 'too far away particles' +! ! also used to determine number of already released shells +! ! time and dtlast are used together to determine 'time passed' +! ! +! call inject_particles(& +! time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,npart,npartoftype,dtinject& +! ) +!end subroutine +! +!subroutine amuse_update_injected() +! use partinject, only:update_injected_particles +! use part, only:npart +! use timestep_ind, only:istepfrac,nbinmax +! use timestep, only:time,dtmax,dt +! implicit none +! integer :: npart_old +! call update_injected_particles(& +! npart_old,npart,istepfrac,nbinmax,time,dtmax,dt,dtinject& +! ) +!end subroutine From 0f8a5837cdc26c6911e6bfb4f650c910df65e7cc Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Tue, 4 Jun 2024 14:03:14 +0200 Subject: [PATCH 08/15] distinguish accretion radius and effective radius for sinks --- src/utils/libphantom-amuse.F90 | 48 +++++++++++++++++++++++++++++----- 1 file changed, 41 insertions(+), 7 deletions(-) diff --git a/src/utils/libphantom-amuse.F90 b/src/utils/libphantom-amuse.F90 index 2e0914237..4f8aa5c2f 100644 --- a/src/utils/libphantom-amuse.F90 +++ b/src/utils/libphantom-amuse.F90 @@ -1112,12 +1112,12 @@ subroutine amuse_new_dm_particle(i, mass, x, y, z, vx, vy, vz) end subroutine subroutine amuse_new_sink_particle(j, mass, x, y, z, vx, vy, vz, & - radius, h_smooth) + radius, accretion_radius, h_smooth) use io, only:fatal - use part, only:nptmass,maxptmass,xyzmh_ptmass,vxyz_ptmass,ihacc,ihsoft + use part, only:nptmass,maxptmass,xyzmh_ptmass,vxyz_ptmass,ihacc,ihsoft,iReff implicit none integer :: i, j - double precision :: mass, x, y, z, vx, vy, vz, radius, h_smooth + double precision :: mass, x, y, z, vx, vy, vz, radius, accretion_radius, h_smooth double precision :: position(3), velocity(3) nptmass = nptmass + 1 @@ -1130,7 +1130,8 @@ subroutine amuse_new_sink_particle(j, mass, x, y, z, vx, vy, vz, & xyzmh_ptmass(2,i) = y xyzmh_ptmass(3,i) = z xyzmh_ptmass(4,i) = mass - xyzmh_ptmass(ihacc,i) = radius + xyzmh_ptmass(iReff,i) = radius + xyzmh_ptmass(ihacc,i) = accretion_radius xyzmh_ptmass(ihsoft,i) = h_smooth vxyz_ptmass(1,i) = vx vxyz_ptmass(2,i) = vy @@ -1338,24 +1339,40 @@ subroutine amuse_get_state_dm(i, mass, x, y, z, vx, vy, vz) write(*,*) 'getting dm ', i end subroutine amuse_get_state_dm -subroutine amuse_get_state_sink(i, mass, x, y, z, vx, vy, vz, radius) +subroutine amuse_get_state_sink(i, mass, x, y, z, vx, vy, vz, radius, accretion_radius) implicit none integer :: i - double precision :: mass, x, y, z, vx, vy, vz, radius + double precision :: mass, x, y, z, vx, vy, vz, radius, accretion_radius call amuse_get_mass(i, mass) call amuse_get_position(i, x, y, z) call amuse_get_velocity(i, vx, vy, vz) call amuse_get_sink_radius(i, radius) write(*,*) 'getting sink ', i, ': radius is ', radius + call amuse_get_sink_accretion_radius(i, accretion_radius) end subroutine amuse_get_state_sink subroutine amuse_get_sink_radius(j, radius) + implicit none + integer :: j + double precision :: radius + call amuse_get_sink_effective_radius(j, radius) +end subroutine amuse_get_sink_radius + +subroutine amuse_get_sink_effective_radius(j, radius) + use part, only:xyzmh_ptmass, iReff + implicit none + integer :: j + double precision :: radius + radius = xyzmh_ptmass(iReff, -j) +end subroutine amuse_get_sink_effective_radius + +subroutine amuse_get_sink_accretion_radius(j, radius) use part, only:xyzmh_ptmass, ihacc implicit none integer :: j double precision :: radius radius = xyzmh_ptmass(ihacc, -j) -end subroutine amuse_get_sink_radius +end subroutine amuse_get_sink_accretion_radius subroutine amuse_get_sink_temperature(j, temperature) use part, only:xyzmh_ptmass, iTeff @@ -1614,6 +1631,15 @@ subroutine amuse_set_state_sink(i, mass, x, y, z, vx, vy, vz, radius) end subroutine subroutine amuse_set_sink_radius(j, radius) + use part, only:xyzmh_ptmass, ihacc, iReff + implicit none + integer :: j + double precision :: radius + call amuse_set_sink_effective_radius(j, radius) + call amuse_set_sink_accretion_radius(j, radius) +end subroutine + +subroutine amuse_set_sink_accretion_radius(j, radius) use part, only:xyzmh_ptmass, ihacc implicit none integer :: j @@ -1621,6 +1647,14 @@ subroutine amuse_set_sink_radius(j, radius) xyzmh_ptmass(ihacc, -j) = radius end subroutine +subroutine amuse_set_sink_effective_radius(j, radius) + use part, only:xyzmh_ptmass, iReff + implicit none + integer :: j + double precision :: radius + xyzmh_ptmass(iReff, -j) = radius +end subroutine + subroutine amuse_set_position(i, x, y, z) use part, only:xyzh,xyzmh_ptmass implicit none From 7c5a63d32f4b7b1adab420a460ff8f03dac02de2 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Mon, 10 Jun 2024 14:47:33 +0200 Subject: [PATCH 09/15] clarification comment --- src/utils/libphantom-amuse.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/utils/libphantom-amuse.F90 b/src/utils/libphantom-amuse.F90 index 4f8aa5c2f..a1054aff4 100644 --- a/src/utils/libphantom-amuse.F90 +++ b/src/utils/libphantom-amuse.F90 @@ -1788,6 +1788,7 @@ subroutine amuse_evolve_model(tmax_in) tlast = time write(*,*) "TIMESTEPPING: evolve from ", time, " to ", tmax + ! Allowing for a shortage of 1% of dtmax to account for floating point differences timestepping: do while (time+0.01 * dtmax < tmax) #ifdef IND_TIMESTEPS From 53c55ff69848ff450a3133578ff0a09d9b2847ea Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Mon, 10 Jun 2024 14:48:10 +0200 Subject: [PATCH 10/15] remove some debugging prints --- src/utils/libphantom-amuse.F90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/utils/libphantom-amuse.F90 b/src/utils/libphantom-amuse.F90 index a1054aff4..c8465787a 100644 --- a/src/utils/libphantom-amuse.F90 +++ b/src/utils/libphantom-amuse.F90 @@ -1347,7 +1347,6 @@ subroutine amuse_get_state_sink(i, mass, x, y, z, vx, vy, vz, radius, accretion_ call amuse_get_position(i, x, y, z) call amuse_get_velocity(i, vx, vy, vz) call amuse_get_sink_radius(i, radius) - write(*,*) 'getting sink ', i, ': radius is ', radius call amuse_get_sink_accretion_radius(i, accretion_radius) end subroutine amuse_get_state_sink @@ -1464,10 +1463,8 @@ subroutine amuse_get_radius(i, radius) integer :: j double precision, intent(out) :: radius if (i == abs(i)) then - write(*,*) "not a sink" call amuse_get_smoothing_length(i, radius) else - write(*,*) "a sink" j = -i call amuse_get_sink_radius(i, radius) endif From d90d9a38e1b090a9cc24d97742d092e574965925 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Fri, 5 Jul 2024 10:36:19 +0200 Subject: [PATCH 11/15] worse indentation, but easier to view diff --- src/main/initial.F90 | 64 ++++++++++++++++++++++---------------------- 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/src/main/initial.F90 b/src/main/initial.F90 index d27f01159..85bf0613a 100644 --- a/src/main/initial.F90 +++ b/src/main/initial.F90 @@ -570,17 +570,17 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread) call init_inject(ierr) if (ierr /= 0) call fatal('initial','error initialising particle injection') if (write_files) then - !rename wind profile filename - inquire(file='wind_profile1D.dat',exist=iexist) - if (iexist) then - i = len(trim(dumpfile)) - if (dumpfile(i-2:i) == 'tmp') then - file1D = dumpfile(1:i-9) // '1D.dat' - else - file1D = dumpfile(1:i-5) // '1D.dat' - endif - call rename('wind_profile1D.dat',trim(file1D)) + !rename wind profile filename + inquire(file='wind_profile1D.dat',exist=iexist) + if (iexist) then + i = len(trim(dumpfile)) + if (dumpfile(i-2:i) == 'tmp') then + file1D = dumpfile(1:i-9) // '1D.dat' + else + file1D = dumpfile(1:i-5) // '1D.dat' endif + call rename('wind_profile1D.dat',trim(file1D)) + endif endif npart_old = npart call inject_particles(time,0.,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& @@ -778,32 +778,32 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread) ! if the input file ends in .tmp or .init ! if (write_files) then - iposinit = index(dumpfile,'.init') - ipostmp = index(dumpfile,'.tmp') - if (iposinit > 0 .or. ipostmp > 0) then + iposinit = index(dumpfile,'.init') + ipostmp = index(dumpfile,'.tmp') + if (iposinit > 0 .or. ipostmp > 0) then #ifdef HDF5 - dumpfileold = trim(dumpfile)//'.h5' + dumpfileold = trim(dumpfile)//'.h5' #else - dumpfileold = dumpfile + dumpfileold = dumpfile #endif - if (iposinit > 0) then - dumpfile = trim(dumpfile(1:iposinit-1)) - else - dumpfile = trim(dumpfile(1:ipostmp-1)) - endif - call write_fulldump(time,trim(dumpfile)) - if (id==master) call write_infile(infile,logfile,evfile,trim(dumpfile),iwritein,iprint) - ! - ! delete temporary dump file - ! - call barrier_mpi() ! Ensure all procs have read temp file before deleting - inquire(file=trim(dumpfileold),exist=iexist) - if (id==master .and. iexist) then - write(iprint,"(/,a,/)") ' ---> DELETING temporary dump file '//trim(dumpfileold)//' <---' - open(unit=idisk1,file=trim(dumpfileold),status='old') - close(unit=idisk1,status='delete') - endif + if (iposinit > 0) then + dumpfile = trim(dumpfile(1:iposinit-1)) + else + dumpfile = trim(dumpfile(1:ipostmp-1)) endif + call write_fulldump(time,trim(dumpfile)) + if (id==master) call write_infile(infile,logfile,evfile,trim(dumpfile),iwritein,iprint) + ! + ! delete temporary dump file + ! + call barrier_mpi() ! Ensure all procs have read temp file before deleting + inquire(file=trim(dumpfileold),exist=iexist) + if (id==master .and. iexist) then + write(iprint,"(/,a,/)") ' ---> DELETING temporary dump file '//trim(dumpfileold)//' <---' + open(unit=idisk1,file=trim(dumpfileold),status='old') + close(unit=idisk1,status='delete') + endif + endif endif ! (write_files) if (id==master) then From f00c75e3e2a37992adacb65ec112a3f2a6be7a03 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Fri, 4 Oct 2024 12:15:41 +0200 Subject: [PATCH 12/15] Modifications to evolve - option not to write files - initialise only on the first call to evol --- src/main/evolve.F90 | 17 +++++++++++++---- src/main/options.f90 | 2 +- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/src/main/evolve.F90 b/src/main/evolve.F90 index 92c22f776..d7e1a9b9f 100644 --- a/src/main/evolve.F90 +++ b/src/main/evolve.F90 @@ -27,6 +27,7 @@ module evolve public :: evol private + logical :: initialized=.false. contains @@ -37,7 +38,7 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag) dtmax_ifactor,dtmax_ifactorWT,dtmax_dratio,check_dtmax_for_decrease,& idtmax_n,idtmax_frac,idtmax_n_next,idtmax_frac_next 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 + 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 @@ -138,6 +139,7 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag) integer :: j,nskip,nskipped,nevwrite_threshold,nskipped_sink,nsinkwrite_threshold character(len=120) :: dumpfile_orig + if (.not. initialized) then tprint = 0. nsteps = 0 nsteplast = 0 @@ -215,6 +217,8 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag) call flush(iprint) + initialised = .true. + endif ! Initialising done ! ! --------------------- main loop ---------------------------------------- ! @@ -387,7 +391,12 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag) if (nskipped >= nevwrite_threshold .or. at_dump_time .or. dt_changed .or. iverbose==5) then nskipped = 0 call get_timings(t1,tcpu1) - call write_evfile(time,dt) + ! 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') @@ -420,7 +429,7 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag) nskipped_sink = nskipped_sink + nskip if (nskipped_sink >= nsinkwrite_threshold .or. at_dump_time .or. dt_changed) then nskipped_sink = 0 - call pt_write_sinkev(nptmass,time,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,fxyz_ptmass_sinksink) + 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 @@ -428,7 +437,7 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag) ! !--write to data file if time is right ! - if (at_dump_time) then + if (at_dump_time .and. write_files) then #ifndef IND_TIMESTEPS ! !--Global timesteps: Decrease dtmax if requested (done in step for individual timesteps) diff --git a/src/main/options.f90 b/src/main/options.f90 index 4041a28ac..f70f43e3c 100644 --- a/src/main/options.f90 +++ b/src/main/options.f90 @@ -172,7 +172,7 @@ subroutine set_default_options ! variable composition use_var_comp = .false. - ! for use as a library + ! enable/disable writing output files write_files = .true. end subroutine set_default_options From 8010c2f9c32925a69db83efc71082a44e537eb64 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Fri, 4 Oct 2024 12:22:20 +0200 Subject: [PATCH 13/15] updates to libphantom-amuse - create, use and update a map to the original particle ids to ensure AMUSE identifies the right particles --- src/utils/libphantom-amuse.F90 | 1232 +++++++++++++++++--------------- 1 file changed, 655 insertions(+), 577 deletions(-) diff --git a/src/utils/libphantom-amuse.F90 b/src/utils/libphantom-amuse.F90 index c8465787a..be8f8167a 100644 --- a/src/utils/libphantom-amuse.F90 +++ b/src/utils/libphantom-amuse.F90 @@ -1,89 +1,54 @@ ! AMUSE interface library for Phantom ! (c) 2019 - 2024 Steven Rieder -! -! Initialize Phantom and set default parameters -! - -subroutine amuse_initialize_code_new() - use dim, only:tagline - use dim, only: maxp_hard - use memory, only:allocate_memory - use mpiutils, only:init_mpi,finalise_mpi - use initial, only:initialise,finalise,startrun,endrun - use io, only:id,master,nprocs,set_io_unit_numbers,die - use evolve, only:evol - use units, only:set_units,utime,umass,udist,unit_density - use physcon, only:gram,seconds,solarm,pc,au - implicit none - integer :: nargs - character(len=120) :: infile,logfile,evfile,dumpfile +module AmusePhantom - id = 0 + use iso_fortran_env + integer(kind=8), allocatable :: amuse_id_lookup(:) + integer :: new_particles_since_last_update = 0 + integer :: particles_added_by_amuse = 0 +contains - call init_mpi(id,nprocs) - call allocate_memory(int(maxp_hard,kind=8)) - call set_io_unit_numbers() - call initialise() - call amuse_set_defaults() ! replaces reading infile - call set_units(dist=1 * au,mass=1.*solarm,G=1.) - -end subroutine amuse_initialize_code_new +subroutine construct_id_lookup() + ! amuse_id_lookup needs to be rebuilt if/when particles are deleted/added/reshuffled + ! easier (though a bit slower) is to do it at the end of an evolve step, and at recommit_particles + use dim, only: maxp, maxp_hard + use part, only: iorig, norig + implicit none + integer :: i, j + write(*,*) "Rebuilding lookup table" + do i=1, maxp + amuse_id_lookup(i) = 0 + enddo + do i=1, norig + j = iorig(i) + if (j > 0) then + amuse_id_lookup(j) = i + endif + !write(*,*) "j = ", j, "i = ", i + !flush(output_unit) + 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() - use dim, only:maxp,maxp_hard,maxvxyzu - use io, only:id,nprocs,iverbose - use mpiutils, only:init_mpi - use memory, only:allocate_memory - use units, only:set_units,utime,umass,udist,unit_density - use physcon, only:gram,seconds,solarm,pc,au - use timestep, only:dtmax,dtextforce - use initial, only:initialise -#ifdef IND_TIMESTEPS - use timestep_ind, only:istepfrac,ibinnow - use part, only:ibin,ibin_old,ibin_wake -#else - use timestep, only:dtcourant,dtforce -#endif + use dim, only:maxp_hard + use allocutils, only:allocate_array + use mpiutils, only:init_mpi + use io, only:id,nprocs implicit none - iverbose=5 - call init_mpi(id,nprocs) - call allocate_memory(int(maxp_hard,kind=8)) - call initialise() - call amuse_set_defaults() - call amuse_set_polyk(0.) - !print*, "maxvxyzu: ", maxvxyzu - !maxvxyzu = 4 - !call set_units(dist=50.*pc,mass=4600.*solarm,G=1.) - !call set_units(dist=1.d20,mass=1.d40,G=1.) - !call set_units(dist=1.*pc,mass=1.*solarm,G=1.) - ! call set_units(dist=1.,mass=1.,time=1.) - - !umass = 1.98892d33 * gram ! AMUSE MSun - !utime = 60 * 60 * 24 * 365.25 * 1d6 * seconds ! 10 Julian Kyr - !call set_units(time=utime,mass=umass,G=1.) - !call set_units(dist=0.1*pc,mass=1.*solarm,G=1.) - call set_units(dist=1 * au,mass=1.*solarm,G=1.) - !print*, "utime/mass/dist/udens: ", utime, umass, udist, unit_density - !dtmax = 0.01 * utime - dtextforce = 1.e-6 -#ifdef IND_TIMESTEPS - ibin(:) = 0 - ibin_old(:) = 0 - ibin_wake(:) = 0 - istepfrac = 0 - ibinnow = 0 -#else - dtcourant = huge(dtcourant) - dtforce = huge(dtforce) -#endif + id = 0 - !call amuse_initialize_wind() + call init_mpi(id,nprocs) + call allocate_array('amuse_id_lookup', amuse_id_lookup, maxp_hard) end subroutine amuse_initialize_code subroutine amuse_set_phantom_option(name, valstring, imatch) ! This subroutine is meant to be a replacement for read_infile + ! It should be kept up to date and support all options listed there, or raise an error if something is not caught. use inject, only:read_options_inject use dust, only:read_options_dust use growth, only:read_options_growth @@ -102,7 +67,7 @@ subroutine amuse_set_phantom_option(name, valstring, imatch) use options, only:alpha,alphaB,alphamax,alphau,twallmax,tolh,rkill,rhofinal_cgs,& psidecayfac,overcleanfac,nmaxdumps,nfulldump,limit_radiation_flux,& ishock_heating,iresistive_heating,ireconav,ipdv_heating,iopacity_type,& - implicit_radiation,hdivbbmax_max,exchange_radiation_energy,calc_erot,& + implicit_radiation,exchange_radiation_energy,calc_erot,& beta,avdecayconst use radiation_implicit, only:tol_rad,itsmax_rad,cv_type use radiation_utils, only:kappa_cgs @@ -119,139 +84,139 @@ subroutine amuse_set_phantom_option(name, valstring, imatch) imatch = .true. select case(trim(name)) case('logfile') - ! ignored + ! ignored case('dumpfile') - ! ignored + ! ignored case('tmax') - read(valstring,*,iostat=ierr) tmax + write(*,*) "Set tmax to ", valstring + read(valstring,*,iostat=ierr) tmax case('dtmax') - read(valstring,*,iostat=ierr) dtmax + read(valstring,*,iostat=ierr) dtmax case('nmax') - read(valstring,*,iostat=ierr) nmax + read(valstring,*,iostat=ierr) nmax case('nout') - read(valstring,*,iostat=ierr) nout + read(valstring,*,iostat=ierr) nout case('nmaxdumps') - read(valstring,*,iostat=ierr) nmaxdumps + read(valstring,*,iostat=ierr) nmaxdumps case('twallmax') - read(valstring,*,iostat=ierr) twallmax + read(valstring,*,iostat=ierr) twallmax case('dtwallmax') - read(valstring,*,iostat=ierr) dtwallmax + read(valstring,*,iostat=ierr) dtwallmax case('iverbose') - read(valstring,*,iostat=ierr) iverbose + read(valstring,*,iostat=ierr) iverbose + write(*,*) "IVERBOSE: ", iverbose case('rhofinal_cgs') - read(valstring,*,iostat=ierr) rhofinal_cgs - incl_runtime2 = .true. + read(valstring,*,iostat=ierr) rhofinal_cgs + incl_runtime2 = .true. case('calc_erot') - read(valstring,*,iostat=ierr) calc_erot - incl_runtime2 = .true. + read(valstring,*,iostat=ierr) calc_erot + incl_runtime2 = .true. case('dtmax_dratio') - read(valstring,*,iostat=ierr) dtmax_dratio - incl_runtime2 = .true. + read(valstring,*,iostat=ierr) dtmax_dratio + incl_runtime2 = .true. case('dtmax_max') - read(valstring,*,iostat=ierr) dtmax_max - if (dtmax_max <= 0.0) dtmax_max = dtmax - ! to prevent comparison errors from round-off - ratio = dtmax_max/dtmax - ratio = int(ratio+0.5)+0.0001 - dtmax_max = dtmax*ratio + read(valstring,*,iostat=ierr) dtmax_max + if (dtmax_max <= 0.0) dtmax_max = dtmax + ! to prevent comparison errors from round-off + ratio = dtmax_max/dtmax + ratio = int(ratio+0.5)+0.0001 + dtmax_max = dtmax*ratio case('dtmax_min') - read(valstring,*,iostat=ierr) dtmax_min - ! to prevent comparison errors from round-off - if (dtmax_min > epsilon(dtmax_min)) then - ratio = dtmax/dtmax_min - ratio = int(ratio+0.5)+0.0001 - dtmax_min = dtmax/ratio - endif + read(valstring,*,iostat=ierr) dtmax_min + ! to prevent comparison errors from round-off + if (dtmax_min > epsilon(dtmax_min)) then + ratio = dtmax/dtmax_min + ratio = int(ratio+0.5)+0.0001 + dtmax_min = dtmax/ratio + endif case('C_cour') - read(valstring,*,iostat=ierr) C_cour + read(valstring,*,iostat=ierr) C_cour case('C_force') - read(valstring,*,iostat=ierr) C_force + read(valstring,*,iostat=ierr) C_force case('tolv') - read(valstring,*,iostat=ierr) tolv + read(valstring,*,iostat=ierr) tolv case('C_ent') - read(valstring,*,iostat=ierr) C_ent + read(valstring,*,iostat=ierr) C_ent case('xtol') - read(valstring,*,iostat=ierr) xtol + read(valstring,*,iostat=ierr) xtol case('ptol') - read(valstring,*,iostat=ierr) ptol + read(valstring,*,iostat=ierr) ptol case('hfact') - read(valstring,*,iostat=ierr) hfact + read(valstring,*,iostat=ierr) hfact case('tolh') - read(valstring,*,iostat=ierr) tolh + read(valstring,*,iostat=ierr) tolh case('rkill') - read(valstring,*,iostat=ierr) rkill + read(valstring,*,iostat=ierr) rkill case('nfulldump') - read(valstring,*,iostat=ierr) nfulldump + read(valstring,*,iostat=ierr) nfulldump case('alpha') - read(valstring,*,iostat=ierr) alpha + read(valstring,*,iostat=ierr) alpha case('alphamax') - read(valstring,*,iostat=ierr) alphamax + read(valstring,*,iostat=ierr) alphamax case('alphau') - read(valstring,*,iostat=ierr) alphau + read(valstring,*,iostat=ierr) alphau case('alphaB') - read(valstring,*,iostat=ierr) alphaB + read(valstring,*,iostat=ierr) alphaB case('psidecayfac') - read(valstring,*,iostat=ierr) psidecayfac + read(valstring,*,iostat=ierr) psidecayfac case('overcleanfac') - read(valstring,*,iostat=ierr) overcleanfac - case('hdivbbmax_max') - read(valstring,*,iostat=ierr) hdivbbmax_max + read(valstring,*,iostat=ierr) overcleanfac case('beta') - read(valstring,*,iostat=ierr) beta + read(valstring,*,iostat=ierr) beta case('ireconav') - read(valstring,*,iostat=ierr) ireconav + read(valstring,*,iostat=ierr) ireconav case('avdecayconst') - read(valstring,*,iostat=ierr) avdecayconst + read(valstring,*,iostat=ierr) avdecayconst case('ipdv_heating') - read(valstring,*,iostat=ierr) ipdv_heating + read(valstring,*,iostat=ierr) ipdv_heating case('ishock_heating') - read(valstring,*,iostat=ierr) ishock_heating + read(valstring,*,iostat=ierr) ishock_heating case('iresistive_heating') - read(valstring,*,iostat=ierr) iresistive_heating + read(valstring,*,iostat=ierr) iresistive_heating case('ien_type') - read(valstring,*,iostat=ierr) ien_type + read(valstring,*,iostat=ierr) ien_type case('irealvisc') - read(valstring,*,iostat=ierr) irealvisc + read(valstring,*,iostat=ierr) irealvisc case('shearparam') - read(valstring,*,iostat=ierr) shearparam + read(valstring,*,iostat=ierr) shearparam case('bulkvisc') - read(valstring,*,iostat=ierr) bulkvisc + read(valstring,*,iostat=ierr) bulkvisc #ifdef MCFOST case('use_mcfost') - read(valstring,*,iostat=ierr) use_mcfost + read(valstring,*,iostat=ierr) use_mcfost case('Voronoi_limits_file') - read(valstring,*,iostat=ierr) Voronoi_limits_file - use_Voronoi_limits_file = .true. + read(valstring,*,iostat=ierr) Voronoi_limits_file + use_Voronoi_limits_file = .true. case('use_mcfost_stars') - read(valstring,*,iostat=ierr) use_mcfost_stellar_parameters + read(valstring,*,iostat=ierr) use_mcfost_stellar_parameters case('mcfost_computes_Lacc') - read(valstring,*,iostat=ierr) mcfost_computes_Lacc + read(valstring,*,iostat=ierr) mcfost_computes_Lacc case('mcfost_uses_PdV') - read(valstring,*,iostat=ierr) mcfost_uses_PdV + read(valstring,*,iostat=ierr) mcfost_uses_PdV case('mcfost_keep_part') - read(valstring,*,iostat=ierr) mcfost_keep_part + read(valstring,*,iostat=ierr) mcfost_keep_part case('ISM') - read(valstring,*,iostat=ierr) ISM + read(valstring,*,iostat=ierr) ISM case('mcfost_dust_subl') - read(valstring,*,iostat=ierr) mcfost_dust_subl + read(valstring,*,iostat=ierr) mcfost_dust_subl #endif case('implicit_radiation') - read(valstring,*,iostat=ierr) implicit_radiation - if (implicit_radiation) store_dust_temperature = .true. + read(valstring,*,iostat=ierr) implicit_radiation + if (implicit_radiation) store_dust_temperature = .true. case('gas-rad_exchange') - read(valstring,*,iostat=ierr) exchange_radiation_energy + read(valstring,*,iostat=ierr) exchange_radiation_energy case('flux_limiter') - read(valstring,*,iostat=ierr) limit_radiation_flux + read(valstring,*,iostat=ierr) limit_radiation_flux case('iopacity_type') - read(valstring,*,iostat=ierr) iopacity_type + read(valstring,*,iostat=ierr) iopacity_type case('cv_type') - read(valstring,*,iostat=ierr) cv_type + read(valstring,*,iostat=ierr) cv_type case('kappa_cgs') - read(valstring,*,iostat=ierr) kappa_cgs + read(valstring,*,iostat=ierr) kappa_cgs case('tol_rad') - read(valstring,*,iostat=ierr) tol_rad + read(valstring,*,iostat=ierr) tol_rad case('itsmax_rad') - read(valstring,*,iostat=ierr) itsmax_rad + read(valstring,*,iostat=ierr) itsmax_rad case default imatch = .false. if (.not.imatch) call read_options_inject(name,valstring,imatch,igotall,ierr) @@ -268,7 +233,7 @@ subroutine amuse_set_phantom_option(name, valstring, imatch) if (.not.imatch) call read_options_gravitationalwaves(name,valstring,imatch,igotall,ierr) if (.not.imatch) call read_options_boundary(name,valstring,imatch,igotall,ierr) end select - if (.not.imatch) write(*,*) "Could not set option ", name + if (.not.imatch) write(*,*) "Could not set option ", name, ", please check if this is a problem!" end subroutine amuse_set_phantom_option subroutine amuse_initialize_wind() @@ -301,39 +266,41 @@ subroutine amuse_initialize_wind() 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 + + call allocate_memory(int(maxp_hard,kind=8)) + call set_io_unit_numbers() + call initialise() + call amuse_set_defaults() ! replaces reading infile + + ! 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 eos, only:ieos,init_eos - !use deriv, only:derivs - !use part, only:npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,Bevol,dBevol,& - ! dustprop,ddustprop,dustfrac,ddustevol,dens,drad,radprop,dustevol,& - ! eos_vars,metrics,pxyzu,rad - !use timestep, only:time,dtmax - !use units, only:udist,utime,umass + use part, only: norig use initial, only:startrun use energies, only: mtot - !use timestep, only:dtmax implicit none - !integer :: ierr - !integer :: nactive - !double precision :: dt, dtnew_first character(len=120) :: infile,logfile,evfile,dumpfile call startrun(infile,logfile,evfile,dumpfile, .true.) - print*, "total mass in code unit: ", mtot - !dtnew_first = dtmax - !dt = 0 - !nactive = npart - print*, "COMMIT_PARTICLES, calling amuse_init_evol" - !call amuse_init_evol() call amuse_evol(.true.) + call construct_id_lookup() + new_particles_since_last_update = norig - particles_added_by_amuse end subroutine amuse_commit_particles subroutine amuse_recommit_particles() use deriv, only:derivs use part, only:npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,Bevol,dBevol,& - dustprop,ddustprop,dustfrac,ddustevol,dens,drad,radprop,dustevol,& + dustprop,filfac,ddustprop,dustfrac,ddustevol,dens,drad,radprop,dustevol,& eos_vars,metrics,pxyzu,rad use timestep, only:time,dtmax implicit none @@ -344,9 +311,9 @@ subroutine amuse_recommit_particles() dt = 0 nactive = npart call derivs(1,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,& - Bevol,dBevol,rad,drad,radprop,dustprop,ddustprop,& - dustevol,ddustevol,dustfrac,eos_vars,time,dt,dtnew_first,pxyzu,dens,metrics) - print*, "RECOMMIT_PARTICLES" + Bevol,dBevol,rad,drad,radprop,dustprop,ddustprop,dustevol,& + ddustevol,filfac,dustfrac,eos_vars,time,dt,dtnew_first,pxyzu,dens,metrics) + call construct_id_lookup() end subroutine amuse_recommit_particles subroutine amuse_cleanup_code() @@ -354,8 +321,16 @@ subroutine amuse_cleanup_code() implicit none call endrun() + if (allocated(amuse_id_lookup)) deallocate(amuse_id_lookup) end subroutine amuse_cleanup_code +subroutine amuse_get_norig(norig_out) + use part, only:norig + implicit none + integer, intent(out) :: norig_out + norig_out = norig +end subroutine amuse_get_norig + ! Get npart subroutine amuse_get_npart(npart_out, nodisabled) use part, only:npart,xyzh @@ -383,9 +358,6 @@ subroutine amuse_set_defaults() implicit none call set_default_options() - ! A few changes to Phantom's defaults - call amuse_set_gamma(1.) - call amuse_set_ieos(1) write_files=.false. end subroutine amuse_set_defaults @@ -398,8 +370,7 @@ subroutine amuse_evol(amuse_initialise) 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,hdivBB_xa,& - compute_energies + 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 @@ -450,10 +421,11 @@ subroutine amuse_evol(amuse_initialise) use io, only:ianalysis #endif use part, only:npart,nptmass,xyzh,vxyzu,fxyzu,fext,divcurlv,massoftype, & - xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,gravity,iboundary, & + 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 + 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 @@ -468,6 +440,7 @@ subroutine amuse_evol(amuse_initialise) #endif implicit none logical, intent(in) :: amuse_initialise + logical :: initialized integer :: flag character(len=256) :: infile @@ -502,12 +475,13 @@ subroutine amuse_evol(amuse_initialise) integer :: j,nskip,nskipped,nevwrite_threshold,nskipped_sink,nsinkwrite_threshold character(len=120) :: dumpfile_orig + initialized = .not. amuse_initialise tzero = time - tlast = time + tlast = time write(*,*) "AMUSE_EVOL called, amuse_initialise=", amuse_initialise write(*,*) "dtlast", dtlast write(*,*) "tzero, tmax: ", tzero, dtmax, tmax - if (amuse_initialise) then + if (.not. initialized) then tprint = 0. nsteps = 0 nsteplast = 0 @@ -532,13 +506,17 @@ subroutine amuse_evol(amuse_initialise) 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(*,*) write(*,*) "***********tlast: ", tlast - write(*,*) dt = dtmax/2.**nbinmax ! use 2.0 here to allow for step too small nmovedtot = 0 tall = 0. @@ -551,18 +529,20 @@ subroutine amuse_evol(amuse_initialise) 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 @@ -582,21 +562,18 @@ subroutine amuse_evol(amuse_initialise) call flush(iprint) -!end subroutine amuse_init_evol - else ! i.e.: .not. amuse_initialise - !tcheck = time - !npart_old = npart - timesubstepping: do while ((time + 0.01 * dtmax < tmax).and.((nsteps < nmax) .or. (nmax < 0)).and.(rhomaxnow*rhofinal1 < 1.0)) - !dtmax = min(dtmax, tmax-time) ! tried this, doesn't work. + 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 - !write(*,*) "INJECTING" 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 @@ -647,7 +624,7 @@ subroutine amuse_evol(amuse_initialise) ! 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,time) + poten,massoftype,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,fxyz_ptmass_sinksink,dptmass,time) endif ! ! Strang splitting: implicit update for half step @@ -738,13 +715,13 @@ subroutine amuse_evol(amuse_initialise) ! !--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 + 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 @@ -758,9 +735,12 @@ subroutine amuse_evol(amuse_initialise) if (nskipped >= nevwrite_threshold .or. at_dump_time .or. dt_changed .or. iverbose==5) then nskipped = 0 call get_timings(t1,tcpu1) - ! We don't want to write the evfile, but we do want to calculate the energies - !call write_evfile(time,dt) - call compute_energies(time) + ! 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') @@ -769,7 +749,7 @@ subroutine amuse_evol(amuse_initialise) call check_conservation_error(mdust(j),mdust_in(j),1.e-1,'dust mass',decrease=.true.) enddo endif - if (mhd) call check_magnetic_stability(hdivBB_xa) + 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)) @@ -793,207 +773,205 @@ subroutine amuse_evol(amuse_initialise) nskipped_sink = nskipped_sink + nskip if (nskipped_sink >= nsinkwrite_threshold .or. at_dump_time .or. dt_changed) then nskipped_sink = 0 - !call pt_write_sinkev(nptmass,time,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,fxyz_ptmass_sinksink) + 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(*,*) - write(*,*) "***********tlast, dtlast: ", tlast, dtlast - write(*,*) + write(*,*) "***********tlast, dtlast, tmax: ", tlast, dtlast, tmax endif ! !--write to data file if time is right ! -! if (at_dump_time) 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 + if (at_dump_time .and. write_files) then +#ifndef IND_TIMESTEPS ! -! !--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 +!--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) ! -! !--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) +!--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) ! -! 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) + call increment_timer(itimer_fromstart,t2-tstart,tcpu2-tcpustart) + + fulldump = (nout <= 0 .and. mod(noutput,nfulldump)==0) .or. (mod(noutput,nout*nfulldump)==0) ! -!#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) +!--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 ! -!#ifndef IND_TIMESTEPS -! !--Implement dynamic boundaries (for global timestepping) -! if (dynamic_bdy) call update_boundaries(nactive,nactive,npart,abortrun_bdy) -!#endif + 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 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 ( (time >= tmax) .or. ( (nmax > 0) .and. (nsteps >= nmax) ) ) fulldump = .true. ! -! 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 +!--flush any buffered warnings to the log file ! -! 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 + if (id==master) call flush_warnings() ! -! twalllast = t2 -! tcpulast = tcpu2 -! do i = 1,ntimers -! call reset_timer(i) -! enddo +!--write dump file ! -! 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 + 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() @@ -1004,16 +982,15 @@ subroutine amuse_evol(amuse_initialise) !--Write out log file prematurely (if requested based upon nstep, walltime) if ( summary_printnow() ) call summary_printout(iprint,nptmass) - enddo timesubstepping + enddo timestepping endif -!end subroutine amuse_new_step 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,npart,npartoftype,xyzh,vxyzu,massoftype + use part, only:igas,norig,npart,npartoftype,xyzh,vxyzu,massoftype use part, only:abundance,iHI,ih2ratio use partinject, only:add_or_update_particle #ifdef IND_TIMESTEPS @@ -1082,7 +1059,7 @@ subroutine amuse_new_sph_particle(i, mass, x, y, z, vx, vy, vz, u, h) print*, x, y, z, vx, vy, vz, mass, ibin(i), twas(i) endif #endif - + particles_added_by_amuse = particles_added_by_amuse + 1 end subroutine amuse_new_sph_particle subroutine amuse_new_dm_particle(i, mass, x, y, z, vx, vy, vz) @@ -1119,12 +1096,16 @@ subroutine amuse_new_sink_particle(j, mass, x, y, z, vx, vy, vz, & 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 with something AMUSE can handle + ! Replace this fatal exception with something AMUSE can handle if (nptmass > maxptmass) call fatal('creating new sink', 'nptmass > maxptmass') i = nptmass + + ! Sink particles are stored in different arrays than other particles. + ! To be able to distinguish the particle indices on the AMUSE side, we give sinks + ! a negative index and other particles a positive index. j = -i + xyzmh_ptmass(:,i) = 0. xyzmh_ptmass(1,i) = x xyzmh_ptmass(2,i) = y @@ -1141,15 +1122,21 @@ 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 :: j + integer :: part_index if (i == abs(i)) then - call kill_particle(i) + write(*,*) "AMUSE killing a gas particle?" + call amuse_get_index(i, part_index) + ! call kill_particle(part_index) else - j = -i + write(*,*) "AMUSE killing a sink particle?" ! Sink particles can't be killed - so we just set its mass to zero - xyzmh_ptmass(4,j) = 0 + ! xyzmh_ptmass(4, -i) = 0 endif + !if (i == abs(i)) then + !else + !endif end subroutine subroutine amuse_get_unit_length(unit_length_out) @@ -1276,14 +1263,42 @@ subroutine amuse_get_time(time_out) time_out = time 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 + dtmax = dt_in + print*, "dtmax: ", dtmax +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 + ! 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! + if (i /= abs(i)) call fatal('get_index', 'get_index is not for sink particles!') + part_index = amuse_id_lookup(i) +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 double precision :: pmassi double precision, intent(out) :: rho - pmassi = massoftype(abs(iphase(i))) - rho = rhoh(xyzh(4,i), pmassi) + call amuse_get_index(i, part_index) + if (part_index == 0) then + rho = 0 + else + pmassi = massoftype(abs(iphase(part_index))) + rho = rhoh(xyzh(4,part_index), pmassi) + endif end subroutine amuse_get_density subroutine amuse_get_pressure(i, p) @@ -1291,36 +1306,46 @@ subroutine amuse_get_pressure(i, p) use eos, only:ieos,equationofstate implicit none integer :: i, eos_type + integer :: part_index double precision :: pmassi, ponrho, rho, spsound, x, y, z double precision, intent(out) :: p real :: tempi - eos_type = ieos - pmassi = massoftype(abs(iphase(i))) - call amuse_get_density(i, rho) - x = xyzh(1,i) - y = xyzh(2,i) - z = xyzh(3,i) - call equationofstate(eos_type,ponrho,spsound,rho,x,y,z,tempi) - p = ponrho * rho + call amuse_get_index(i, part_index) + if (part_index == 0) then + p = 0 + else + eos_type = ieos + pmassi = massoftype(abs(iphase(part_index))) + call amuse_get_density(part_index, rho) + x = xyzh(1,part_index) + y = xyzh(2,part_index) + z = xyzh(3,part_index) + call equationofstate(eos_type,ponrho,spsound,rho,x,y,z,tempi) + p = ponrho * rho + endif end subroutine amuse_get_pressure subroutine amuse_get_mass(i, part_mass) use part, only:iphase,massoftype,xyzmh_ptmass implicit none double precision, intent(out) :: part_mass - integer, intent(in) :: i - integer :: j + integer, intent(inout) :: i + integer :: part_index if (i == abs(i)) then - part_mass = massoftype(abs(iphase(i))) + call amuse_get_index(i, part_index) + if (part_index == 0) then + part_mass = 0 + else + part_mass = massoftype(abs(iphase(part_index))) + endif else - j = -i - part_mass = xyzmh_ptmass(4,j) + part_mass = xyzmh_ptmass(4, -i) endif end subroutine amuse_get_mass subroutine amuse_get_state_gas(i, mass, x, y, z, vx, vy, vz, u, h) implicit none - integer :: i + integer, 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) @@ -1331,8 +1356,8 @@ end subroutine amuse_get_state_gas subroutine amuse_get_state_dm(i, mass, x, y, z, vx, vy, vz) implicit none - integer :: i - double precision :: mass, x, y, z, vx, vy, vz + integer, 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) call amuse_get_velocity(i, vx, vy, vz) @@ -1341,8 +1366,8 @@ 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 :: i - double precision :: mass, x, y, z, vx, vy, vz, radius, accretion_radius + integer, 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) call amuse_get_velocity(i, vx, vy, vz) @@ -1350,96 +1375,154 @@ subroutine amuse_get_state_sink(i, mass, x, y, z, vx, vy, vz, radius, accretion_ call amuse_get_sink_accretion_radius(i, accretion_radius) end subroutine amuse_get_state_sink -subroutine amuse_get_sink_radius(j, radius) +subroutine amuse_set_state_gas(i, mass, x, y, z, vx, vy, vz, u, h) + implicit none + integer :: 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) + call amuse_set_velocity(i, vx, vy, vz) + call amuse_set_internal_energy(i, u) + call amuse_set_smoothing_length(i, h) +end subroutine + +subroutine amuse_set_state_dm(i, mass, x, y, z, vx, vy, vz) + implicit none + integer :: i + double precision :: mass, x, y, z, vx, vy, vz + call amuse_set_mass(i, mass) + call amuse_set_position(i, x, y, z) + call amuse_set_velocity(i, vx, vy, vz) +end subroutine + +subroutine amuse_set_state_sink(i, mass, x, y, z, vx, vy, vz, radius, accretion_radius) implicit none - integer :: j + integer :: 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) + call amuse_set_velocity(i, vx, vy, vz) + call amuse_set_sink_radius(i, radius) + call amuse_set_sink_accretion_radius(i, accretion_radius) +end subroutine + +subroutine amuse_get_sink_radius(i, radius) + implicit none + integer :: i double precision :: radius - call amuse_get_sink_effective_radius(j, radius) + call amuse_get_sink_effective_radius(i, radius) end subroutine amuse_get_sink_radius -subroutine amuse_get_sink_effective_radius(j, radius) +subroutine amuse_set_sink_radius(i, radius) + use part, only:xyzmh_ptmass, ihacc, iReff + implicit none + integer :: i + double precision :: radius + call amuse_set_sink_effective_radius(i, radius) + call amuse_set_sink_accretion_radius(i, radius) +end subroutine + +subroutine amuse_get_sink_effective_radius(i, radius) use part, only:xyzmh_ptmass, iReff implicit none - integer :: j + integer :: i double precision :: radius - radius = xyzmh_ptmass(iReff, -j) + radius = xyzmh_ptmass(iReff, -i) end subroutine amuse_get_sink_effective_radius -subroutine amuse_get_sink_accretion_radius(j, radius) +subroutine amuse_get_sink_accretion_radius(i, radius) use part, only:xyzmh_ptmass, ihacc implicit none - integer :: j + integer :: i double precision :: radius - radius = xyzmh_ptmass(ihacc, -j) + radius = xyzmh_ptmass(ihacc, -i) end subroutine amuse_get_sink_accretion_radius -subroutine amuse_get_sink_temperature(j, temperature) +subroutine amuse_get_sink_temperature(i, temperature) use part, only:xyzmh_ptmass, iTeff implicit none - integer, intent(in) :: j + integer, intent(in) :: i double precision :: temperature - temperature = xyzmh_ptmass(iTeff, -j) + temperature = xyzmh_ptmass(iTeff, -i) end subroutine -subroutine amuse_get_sink_luminosity(j, luminosity) +subroutine amuse_get_sink_luminosity(i, luminosity) use part, only:xyzmh_ptmass, iLum implicit none - integer, intent(in) :: j + integer, intent(in) :: i double precision :: luminosity - luminosity = xyzmh_ptmass(iLum, -j) + luminosity = xyzmh_ptmass(iLum, -i) end subroutine subroutine amuse_get_position(i, x, y, z) use part, only:xyzh,xyzmh_ptmass implicit none - integer, intent(in) :: i - integer :: j + integer, intent(inout) :: i + integer :: part_index double precision, intent(out) :: x, y, z if (i == abs(i)) then - x = xyzh(1, i) - y = xyzh(2, i) - z = xyzh(3, i) + call amuse_get_index(i, part_index) + if (part_index == 0) then + x = 0 + y = 0 + z = 0 + else + x = xyzh(1, part_index) + y = xyzh(2, part_index) + z = xyzh(3, part_index) + endif else - j = -i - x = xyzmh_ptmass(1, j) - y = xyzmh_ptmass(2, j) - z = xyzmh_ptmass(3, j) + x = xyzmh_ptmass(1, -i) + y = xyzmh_ptmass(2, -i) + z = xyzmh_ptmass(3, -i) endif end subroutine amuse_get_position subroutine amuse_get_velocity(i, vx, vy, vz) use part, only:vxyzu,vxyz_ptmass implicit none - integer, intent(in) :: i - integer :: j + integer, intent(inout) :: i + integer :: part_index double precision, intent(out) :: vx, vy, vz if (i == abs(i)) then - vx = vxyzu(1, i) - vy = vxyzu(2, i) - vz = vxyzu(3, i) + call amuse_get_index(i, part_index) + if (part_index == 0) then + vx = 0 + vy = 0 + vz = 0 + else + vx = vxyzu(1, part_index) + vy = vxyzu(2, part_index) + vz = vxyzu(3, part_index) + endif else - j = -i - vx = vxyz_ptmass(1, j) - vy = vxyz_ptmass(2, j) - vz = vxyz_ptmass(3, j) + vx = vxyz_ptmass(1, -i) + vy = vxyz_ptmass(2, -i) + vz = vxyz_ptmass(3, -i) endif end subroutine amuse_get_velocity subroutine amuse_get_acceleration(i, fx, fy, fz) use part, only:fxyzu,fxyz_ptmass implicit none - integer, intent(in) :: i - integer :: j + integer, intent(inout) :: i + integer :: part_index double precision, intent(out) :: fx, fy, fz if (i == abs(i)) then - fx = fxyzu(1, i) - fy = fxyzu(2, i) - fz = fxyzu(3, i) + call amuse_get_index(i, part_index) + if (part_index == 0) then + fx = 0 + fy = 0 + fz = 0 + else + fx = fxyzu(1, part_index) + fy = fxyzu(2, part_index) + fz = fxyzu(3, part_index) + endif else - j = -i - fx = fxyz_ptmass(1, j) - fy = fxyz_ptmass(2, j) - fz = fxyz_ptmass(3, j) + fx = fxyz_ptmass(1, -i) + fy = fxyz_ptmass(2, -i) + fz = fxyz_ptmass(3, -i) endif end subroutine amuse_get_acceleration @@ -1447,25 +1530,27 @@ subroutine amuse_get_smoothing_length(i, h) use part, only:xyzh,xyzmh_ptmass,ihsoft implicit none integer, intent(in) :: i - integer :: j + integer :: part_index double precision, intent(out) :: h if (i == abs(i)) then - h = xyzh(4, i) + call amuse_get_index(i, part_index) + if (part_index == 0) then + h = 0 + else + h = xyzh(4, part_index) + endif else - j = -i - h = xyzmh_ptmass(ihsoft,j) + h = xyzmh_ptmass(ihsoft, -i) endif end subroutine amuse_get_smoothing_length subroutine amuse_get_radius(i, radius) implicit none integer, intent(in) :: i - integer :: j double precision, intent(out) :: radius if (i == abs(i)) then call amuse_get_smoothing_length(i, radius) else - j = -i call amuse_get_sink_radius(i, radius) endif end subroutine amuse_get_radius @@ -1475,10 +1560,13 @@ subroutine amuse_get_internal_energy(i, u) use part, only:vxyzu implicit none integer, intent(in) :: i + integer :: part_index double precision, intent(out) :: u - - if (maxvxyzu >= 4) then - u = vxyzu(4, i) + call amuse_get_index(i, part_index) + if (part_index == 0) then + u = 0 + else if (maxvxyzu >= 4) then + u = vxyzu(4, part_index) else u = 0 endif @@ -1488,99 +1576,125 @@ subroutine amuse_set_hi_abundance(i, hi_abundance) use part, only:abundance,iHI implicit none integer, intent(in) :: i + integer :: part_index double precision, intent(in) :: hi_abundance + call amuse_get_index(i, part_index) - abundance(iHI, i) = hi_abundance + abundance(iHI, part_index) = hi_abundance end subroutine subroutine amuse_get_hi_abundance(i, hi_abundance) use part, only:abundance,iHI implicit none integer, intent(in) :: i + integer :: part_index double precision, intent(out) :: hi_abundance - - hi_abundance = abundance(iHI, i) + call amuse_get_index(i, part_index) + if (part_index == 0) then + hi_abundance = 0 + else + hi_abundance = abundance(iHI, part_index) + endif end subroutine subroutine amuse_set_proton_abundance(i, proton_abundance) use part, only:abundance,iproton implicit none integer, intent(in) :: i + integer :: part_index double precision, intent(in) :: proton_abundance + call amuse_get_index(i, part_index) - abundance(iproton, i) = proton_abundance + abundance(iproton, part_index) = proton_abundance end subroutine subroutine amuse_get_proton_abundance(i, proton_abundance) use part, only:abundance,iproton implicit none integer, intent(in) :: i + integer :: part_index double precision, intent(out) :: proton_abundance - - proton_abundance = abundance(iproton, i) + call amuse_get_index(i, part_index) + if (part_index == 0) then + proton_abundance = 0 + else + proton_abundance = abundance(iproton, part_index) + endif end subroutine subroutine amuse_set_electron_abundance(i, electron_abundance) use part, only:abundance,ielectron implicit none integer, intent(in) :: i + integer :: part_index double precision, intent(in) :: electron_abundance + call amuse_get_index(i, part_index) - abundance(ielectron, i) = electron_abundance + abundance(ielectron, part_index) = electron_abundance end subroutine subroutine amuse_get_electron_abundance(i, electron_abundance) use part, only:abundance,ielectron implicit none integer, intent(in) :: i + integer :: part_index double precision, intent(out) :: electron_abundance - - electron_abundance = abundance(ielectron, i) + call amuse_get_index(i, part_index) + if (part_index == 0) then + electron_abundance = 0 + else + electron_abundance = abundance(ielectron, part_index) + endif end subroutine subroutine amuse_set_co_abundance(i, co_abundance) use part, only:abundance,ico implicit none integer, intent(in) :: i + integer :: part_index double precision, intent(in) :: co_abundance + call amuse_get_index(i, part_index) - abundance(ico, i) = co_abundance + abundance(ico, part_index) = co_abundance end subroutine subroutine amuse_get_co_abundance(i, co_abundance) use part, only:abundance,ico implicit none integer, intent(in) :: i + integer :: part_index double precision, intent(out) :: co_abundance - - co_abundance = abundance(ico, i) + call amuse_get_index(i, part_index) + if (part_index == 0) then + co_abundance = 0 + else + co_abundance = abundance(ico, part_index) + endif end subroutine subroutine amuse_set_h2ratio(i, h2ratio) use part, only:abundance,iHI,ih2ratio implicit none integer, intent(in) :: i + integer :: part_index double precision, intent(in) :: h2ratio + call amuse_get_index(i, part_index) - abundance(ih2ratio, i) = h2ratio + abundance(ih2ratio, part_index) = h2ratio end subroutine subroutine amuse_get_h2ratio(i, h2ratio) use part, only:abundance,iHI,ih2ratio implicit none integer, intent(in) :: i + integer :: part_index double precision, intent(out) :: h2ratio - - h2ratio = abundance(ih2ratio, i) -end subroutine - -subroutine amuse_set_time_step(dt_in) - use units, only:utime - use timestep, only:dtmax - implicit none - double precision, intent(in) :: dt_in - dtmax = dt_in - print*, "dtmax: ", dtmax + call amuse_get_index(i, part_index) + if (part_index == 0) then + h2ratio = 0 + else + h2ratio = abundance(ih2ratio, part_index) + endif end subroutine subroutine amuse_set_mass(i, part_mass) @@ -1588,85 +1702,46 @@ subroutine amuse_set_mass(i, part_mass) implicit none double precision, intent(in) :: part_mass integer, intent(in) :: i - integer :: j + integer :: part_index if (i == abs(i)) then - massoftype(abs(iphase(i))) = part_mass + call amuse_get_index(i, part_index) + massoftype(abs(iphase(part_index))) = part_mass else - j = -i - xyzmh_ptmass(4,j) = part_mass + xyzmh_ptmass(4,part_index) = part_mass endif end subroutine -subroutine amuse_set_state_gas(i, mass, x, y, z, vx, vy, vz, u, h) - implicit none - integer :: 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) - call amuse_set_velocity(i, vx, vy, vz) - call amuse_set_internal_energy(i, u) - call amuse_set_smoothing_length(i, h) -end subroutine - -subroutine amuse_set_state_dm(i, mass, x, y, z, vx, vy, vz) - implicit none - integer :: i - double precision :: mass, x, y, z, vx, vy, vz - call amuse_set_mass(i, mass) - call amuse_set_position(i, x, y, z) - call amuse_set_velocity(i, vx, vy, vz) -end subroutine - -subroutine amuse_set_state_sink(i, mass, x, y, z, vx, vy, vz, radius) - implicit none - integer :: i - double precision :: mass, x, y, z, vx, vy, vz, radius - call amuse_set_mass(i, mass) - call amuse_set_position(i, x, y, z) - call amuse_set_velocity(i, vx, vy, vz) - call amuse_set_sink_radius(i, radius) -end subroutine - -subroutine amuse_set_sink_radius(j, radius) - use part, only:xyzmh_ptmass, ihacc, iReff - implicit none - integer :: j - double precision :: radius - call amuse_set_sink_effective_radius(j, radius) - call amuse_set_sink_accretion_radius(j, radius) -end subroutine - -subroutine amuse_set_sink_accretion_radius(j, radius) +subroutine amuse_set_sink_accretion_radius(i, radius) use part, only:xyzmh_ptmass, ihacc implicit none - integer :: j + integer :: i double precision :: radius - xyzmh_ptmass(ihacc, -j) = radius + xyzmh_ptmass(ihacc, i) = radius end subroutine -subroutine amuse_set_sink_effective_radius(j, radius) +subroutine amuse_set_sink_effective_radius(i, radius) use part, only:xyzmh_ptmass, iReff implicit none - integer :: j + integer :: i double precision :: radius - xyzmh_ptmass(iReff, -j) = radius + xyzmh_ptmass(iReff, i) = radius end subroutine subroutine amuse_set_position(i, x, y, z) use part, only:xyzh,xyzmh_ptmass implicit none integer, intent(in) :: i - integer :: j + integer :: part_index double precision, intent(in) :: x, y, z if (i == abs(i)) then - xyzh(1, i) = x - xyzh(2, i) = y - xyzh(3, i) = z + call amuse_get_index(i, part_index) + xyzh(1, part_index) = x + xyzh(2, part_index) = y + xyzh(3, part_index) = z else - j = -i - xyzmh_ptmass(1, j) = x - xyzmh_ptmass(2, j) = y - xyzmh_ptmass(3, j) = z + xyzmh_ptmass(1, -i) = x + xyzmh_ptmass(2, -i) = y + xyzmh_ptmass(3, -i) = z endif end subroutine @@ -1674,17 +1749,17 @@ subroutine amuse_set_velocity(i, vx, vy, vz) use part, only:vxyzu,vxyz_ptmass implicit none integer, intent(in) :: i - integer :: j + integer :: part_index double precision, intent(in) :: vx, vy, vz if (i == abs(i)) then - vxyzu(1, i) = vx - vxyzu(2, i) = vy - vxyzu(3, i) = vz + call amuse_get_index(i, part_index) + vxyzu(1, part_index) = vx + vxyzu(2, part_index) = vy + vxyzu(3, part_index) = vz else - j = -i - vxyz_ptmass(1, j) = vx - vxyz_ptmass(2, j) = vy - vxyz_ptmass(3, j) = vz + vxyz_ptmass(1, -i) = vx + vxyz_ptmass(2, -i) = vy + vxyz_ptmass(3, -i) = vz endif end subroutine @@ -1692,43 +1767,41 @@ subroutine amuse_set_smoothing_length(i, h) use part, only:xyzh,xyzmh_ptmass,ihsoft implicit none integer, intent(in) :: i - integer :: j + integer :: part_index double precision, intent(in) :: h if (i == abs(i)) then - xyzh(4, i) = h + call amuse_get_index(i, part_index) + xyzh(4, part_index) = h else - j = -i - xyzmh_ptmass(ihsoft, j) = h + xyzmh_ptmass(ihsoft, -i) = h endif end subroutine subroutine amuse_set_radius(i, radius) implicit none integer, intent(in) :: i - integer :: j double precision :: radius if (i == abs(i)) then call amuse_set_smoothing_length(i, radius) else - j = -i - call amuse_set_sink_radius(j, radius) + call amuse_set_sink_radius(i, radius) endif end subroutine -subroutine amuse_set_sink_temperature(j, temperature) +subroutine amuse_set_sink_temperature(i, temperature) use part, only:xyzmh_ptmass, iTeff implicit none - integer, intent(in) :: j + integer, intent(in) :: i double precision :: temperature - xyzmh_ptmass(iTeff, -j) = temperature + xyzmh_ptmass(iTeff, -i) = temperature end subroutine -subroutine amuse_set_sink_luminosity(j, luminosity) +subroutine amuse_set_sink_luminosity(i, luminosity) use part, only:xyzmh_ptmass, iLum implicit none - integer, intent(in) :: j + integer, intent(in) :: i double precision :: luminosity - xyzmh_ptmass(iLum, -j) = luminosity + xyzmh_ptmass(iLum, -i) = luminosity end subroutine subroutine amuse_set_internal_energy(i, u) @@ -1737,9 +1810,11 @@ subroutine amuse_set_internal_energy(i, u) use timestep, only:dtextforce implicit none integer, intent(in) :: i + integer :: part_index double precision, intent(in) :: u + call amuse_get_index(i, part_index) if (maxvxyzu >= 4) then - vxyzu(4, i) = u + vxyzu(4, part_index) = u endif ! Changing temperature -> better use a small cooling step dtextforce = 1.e-8 @@ -1747,7 +1822,6 @@ subroutine amuse_set_internal_energy(i, u) subroutine amuse_evolve_model(tmax_in) use timestep, only:tmax, time, dt, dtmax, rhomaxnow, dtlast - ! use evolvesplit, only:init_step, finalize_step #ifdef IND_TIMESTEPS use timestep_ind, only:istepfrac #endif @@ -1758,11 +1832,13 @@ subroutine amuse_evolve_model(tmax_in) #endif use options, only:rhofinal1 use ptmass, only:rho_crit - use part, only:npart + use part, only:npart, norig use step_lf_global, only:init_step use part, only: xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass implicit none + integer :: number_of_particles_at_start + integer :: number_of_particles_at_finish logical :: amuse_initialise double precision, intent(in) :: tmax_in logical :: maximum_density_reached @@ -1779,6 +1855,7 @@ subroutine amuse_evolve_model(tmax_in) dtinject = huge(dtinject) ! dtlast = 0 nbinmax = 0 + number_of_particles_at_start = norig tmax = tmax_in ! - epsilon(tmax_in) !dtmax = (tmax - time) @@ -1790,15 +1867,14 @@ subroutine amuse_evolve_model(tmax_in) #ifdef IND_TIMESTEPS istepfrac = 0 - !print*, "*****init timestep" - !call init_step(npart, time, dtmax) #endif - !print*, "*****new_step - Time; tmax: ", time, tmax - !call amuse_new_step(tlast) amuse_initialise = .false. call amuse_evol(amuse_initialise) - !print*, "*****new_step done - Time; tmax: ", time, tmax enddo timestepping + + call construct_id_lookup() + number_of_particles_at_finish = norig + new_particles_since_last_update = new_particles_since_last_update + number_of_particles_at_finish - number_of_particles_at_start end subroutine ! @@ -2304,3 +2380,5 @@ subroutine amuse_get_udist(udist_out) ! npart_old,npart,istepfrac,nbinmax,time,dtmax,dt,dtinject& ! ) !end subroutine + +end module From 304f743d5d133c008cd673763039f7f681ce1dfb Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Fri, 4 Oct 2024 13:38:30 +0200 Subject: [PATCH 14/15] corrections --- src/main/evolve.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/main/evolve.F90 b/src/main/evolve.F90 index d7e1a9b9f..b7a3801f7 100644 --- a/src/main/evolve.F90 +++ b/src/main/evolve.F90 @@ -44,7 +44,7 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag) check_magnetic_stability use dim, only:maxvxyzu,mhd,periodic,idumpfile use fileutils, only:getnextfilename - use options, only:nfulldump,twallmax,nmaxdumps,rhofinal1,iexternalforce,rkill + use options, only:nfulldump,twallmax,nmaxdumps,rhofinal1,iexternalforce,rkill,write_files use readwrite_infile, only:write_infile use readwrite_dumps, only:write_smalldump,write_fulldump use step_lf_global, only:step @@ -217,7 +217,7 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag) call flush(iprint) - initialised = .true. + initialized = .true. endif ! Initialising done ! ! --------------------- main loop ---------------------------------------- From d39f1d6cc8218800187ccd104573ce26142f898e Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Fri, 18 Oct 2024 12:17:01 +0200 Subject: [PATCH 15/15] 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