Libphantom-amuse

build/Makefile

@@ -859,6 +859,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
 

build/Makefile_setups

@@ -1122,3 +1122,20 @@ ifeq ($(SETUP), radiativebox)
     RADIATION=yes
     PERIODIC=yes
 endif
+
+ifeq ($(SETUP), amuse)
+    FPPFLAGS=-DAMUSE
+    # H2CHEM=yes
+    ISOTHERMAL=no
+    GRAVITY=yes
+    # IND_TIMESTEPS=yes
+    IND_TIMESTEPS=no
+    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

src/main/evolve.F90

@@ -27,6 +27,7 @@ module evolve
  public :: evol
 
  private
+ logical :: initialized=.false.
 
 contains
 
@@ -37,13 +38,13 @@ 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
  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
@@ -138,10 +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)
@@ -215,6 +217,8 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag)
 
  call flush(iprint)
 
+ initialized = .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,15 +429,15 @@ 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
     endif
 !
 !--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)

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

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, &
@@ -569,6 +569,7 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
 #ifdef INJECT_PARTICLES
  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
@@ -580,6 +581,7 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
     endif
     call rename('wind_profile1D.dat',trim(file1D))
  endif
+ endif
  npart_old = npart
  call inject_particles(time,0.,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,&
                        npart,npart_old,npartoftype,dtinject)
@@ -656,22 +658,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
@@ -775,6 +777,7 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
 !--write initial conditions to output file
 !  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
@@ -801,6 +804,7 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread)
        close(unit=idisk1,status='delete')
     endif
  endif
+ endif ! (write_files)
 
  if (id==master) then
     call flush_warnings()

src/main/options.f90

@@ -55,10 +55,13 @@ module options
  real(kind=sp), 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
@@ -169,6 +172,9 @@ subroutine set_default_options
  ! variable composition
  use_var_comp = .false.
 
+ ! enable/disable writing output files
+ write_files = .true.
+
 end subroutine set_default_options
 
 end module options

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
@@ -1706,6 +1707,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
@@ -1782,6 +1784,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

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
 

src/main/wind.F90

@@ -929,6 +929,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
@@ -953,10 +954,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
@@ -987,7 +989,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
@@ -1007,7 +1009,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)