Merge pull request #544 from danieljprice/grstar-disc

(#543) MAXP should not be hardwired in the SETUP block

build/Makefile_setups

@@ -455,7 +455,6 @@ ifeq ($(SETUP), nimhdshock)
     NONIDEALMHD=yes
     KERNEL=WendlandC4
     ISOTHERMAL=yes
-    MAXP=6000000
     KNOWN_SETUP=yes
 endif
 
@@ -738,7 +737,6 @@ ifeq ($(SETUP), star)
     MODFILE= utils_binary.f90 set_binary.f90 moddump_binary.f90
     ANALYSIS= ${SRCNIMHD} utils_summary.o utils_omp.o ptmass.o energies.o analysis_common_envelope.f90
     KNOWN_SETUP=yes
-    MAXP=10000000
     GRAVITY=yes
 endif
 
@@ -751,7 +749,6 @@ ifeq ($(SETUP), grstar)
     MODFILE= moddump_tidal.f90
     ANALYSIS= ${SRCNIMHD} utils_summary.o utils_omp.o ptmass.o energies.o analysis_common_envelope.f90
     KNOWN_SETUP=yes
-    MAXP=100000000
     GRAVITY=yes
 endif
 
@@ -761,7 +758,6 @@ ifeq ($(SETUP), radstar)
     MODFILE= utils_binary.f90 set_binary.f90 moddump_binary.f90
     ANALYSIS= ${SRCNIMHD} utils_summary.o utils_omp.o ptmass.o energies.o analysis_common_envelope.f90
     KNOWN_SETUP=yes
-    MAXP=10000000
     GRAVITY=yes
     RADIATION=yes
 endif
@@ -773,7 +769,6 @@ ifeq ($(SETUP), dustystar)
     MODFILE= utils_binary.f90 set_binary.f90 moddump_binary.f90
     ANALYSIS= ${SRCNIMHD} utils_summary.o utils_omp.o ptmass.o energies.o analysis_common_envelope.f90 dust_formation.f90
     KNOWN_SETUP=yes
-    MAXP=10000000
     GRAVITY=yes
     SINK_RADIATION=yes
 endif

scripts/buildbot.sh

@@ -184,9 +184,10 @@ check_phantomsetup ()
    #
    myinput="\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n";
    prefix="myrun";
+   flags="--np=1000"
    echo -e "$myinput" > myinput.txt;
    sed '/-e/d' myinput.txt > mycleanin.txt
-   ./phantomsetup $prefix < mycleanin.txt > /dev/null; err=$?;
+   ./phantomsetup $prefix $flags < mycleanin.txt > /dev/null; err=$?;
    if [ $err -eq 0 ]; then
       print_result "runs" $pass;
    else
@@ -197,8 +198,8 @@ check_phantomsetup ()
    # run phantomsetup up to 3 times to successfully create/rewrite the .setup file
    #
    infile="${prefix}.in"
-   ./phantomsetup $prefix < myinput.txt > /dev/null;
-   ./phantomsetup $prefix < myinput.txt > /dev/null;
+   ./phantomsetup $prefix $flags < myinput.txt > /dev/null;
+   ./phantomsetup $prefix $flags < myinput.txt > /dev/null;
    if [ -e "$prefix.setup" ]; then
       print_result "creates .setup file" $pass;
       #test_setupfile_options "$prefix" "$prefix.setup" $infile;

src/main/checksetup.f90

@@ -1020,13 +1020,17 @@ subroutine check_setup_radiation(npart,nerror,nwarn,radprop,rad)
 end subroutine check_setup_radiation
 
 subroutine check_vdep_extf(nwarn,iexternalforce)
- use externalforces, only: is_velocity_dependent
- use ptmass, only : use_fourthorder
+ use externalforces, only:is_velocity_dependent
+ use ptmass,         only:use_fourthorder
+ use dim,            only:gr
  integer, intent(inout) :: nwarn
  integer, intent(in)    :: iexternalforce
- if (is_velocity_dependent(iexternalforce) .and. use_fourthorder) then
-    print "(/,a,/)","Warning: velocity dependant external forces are not compatible with FSI switch back to Leapfrog..."
-    nwarn = nwarn + 1
+
+ if (iexternalforce > 0 .and. is_velocity_dependent(iexternalforce) .and. use_fourthorder) then
+    if (.not.gr) then ! do not give the warning in GR, just do it...
+       print "(/,1x,a,/)"," Warning: Switching to Leapfrog integrator for velocity-dependent external forces..."
+       nwarn = nwarn + 1
+    endif
     use_fourthorder = .false.
  endif
 

src/main/deriv.F90

@@ -58,7 +58,7 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,&
  use forces,         only:force
  use part,           only:mhd,gradh,alphaind,igas,iradxi,ifluxx,ifluxy,ifluxz,ithick
  use derivutils,     only:do_timing
- use cons2prim,      only:cons2primall,cons2prim_everything,prim2consall
+ use cons2prim,      only:cons2primall,cons2prim_everything
  use metric_tools,   only:init_metric
  use radiation_implicit, only:do_radiation_implicit,ierr_failed_to_converge
  use options,        only:implicit_radiation,implicit_radiation_store_drad,use_porosity
@@ -115,7 +115,6 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,&
     if (gr) then
        ! Recalculate the metric after moving particles to their new tasks
        call init_metric(npart,xyzh,metrics)
-       !call prim2consall(npart,xyzh,metrics,vxyzu,dens,pxyzu,use_dens=.false.)
     endif
 
     if (nptmass > 0 .and. periodic) call ptmass_boundary_crossing(nptmass,xyzmh_ptmass)
@@ -221,11 +220,13 @@ end subroutine derivs
 !+
 !--------------------------------------
 subroutine get_derivs_global(tused,dt_new,dt)
- use part,   only:npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,&
-                Bevol,dBevol,rad,drad,radprop,dustprop,ddustprop,filfac,&
-                dustfrac,ddustevol,eos_vars,pxyzu,dens,metrics,dustevol
- use timing, only:printused,getused
- use io,     only:id,master
+ use part,         only:npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,&
+                        Bevol,dBevol,rad,drad,radprop,dustprop,ddustprop,filfac,&
+                        dustfrac,ddustevol,eos_vars,pxyzu,dens,metrics,dustevol,gr
+ use timing,       only:printused,getused
+ use io,           only:id,master
+ use cons2prim,    only:prim2consall
+ use metric_tools, only:init_metric
  real(kind=4), intent(out), optional :: tused
  real,         intent(out), optional :: dt_new
  real,         intent(in), optional  :: dt  ! optional argument needed to test implicit radiation routine
@@ -237,6 +238,13 @@ subroutine get_derivs_global(tused,dt_new,dt)
  dti = 0.
  if (present(dt)) dti = dt
  call getused(t1)
+ ! update conserved quantities in the GR code
+ if (gr) then
+    call init_metric(npart,xyzh,metrics)
+    call prim2consall(npart,xyzh,metrics,vxyzu,dens,pxyzu,use_dens=.false.)
+ endif
+
+ ! evaluate derivatives
  call derivs(1,npart,npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,Bevol,dBevol,&
              rad,drad,radprop,dustprop,ddustprop,dustevol,ddustevol,filfac,dustfrac,&
              eos_vars,time,dti,dtnew,pxyzu,dens,metrics)

src/main/utils_infiles.f90

@@ -443,11 +443,17 @@ subroutine read_inopt_int(ival,tag,db,err,errcount,min,max)
  if (ierr==0) then
     if (present(min)) then
        write(chmin,"(g10.0)") min
-       if (ival < min) ierr = ierr_rangemin
+       if (ival < min) then
+          ierr = ierr_rangemin
+          ival = min
+       endif
     endif
     if (present(max)) then
        write(chmax,"(g10.0)") max
-       if (ival > max) ierr = ierr_rangemax
+       if (ival > max) then
+          ierr = ierr_rangemax
+          ival = max
+       endif
     endif
  endif
 
@@ -493,11 +499,17 @@ subroutine read_inopt_real(val,tag,db,err,errcount,min,max)
  if (ierr==0) then
     if (present(min)) then
        write(chmin,"(g13.4)") min
-       if (val < min) ierr = ierr_rangemin
+       if (val < min) then
+          ierr = ierr_rangemin
+          val = min
+       endif
     endif
     if (present(max)) then
        write(chmax,"(g13.4)") max
-       if (val > max) ierr = ierr_rangemax
+       if (val > max) then
+          ierr = ierr_rangemax
+          val = max
+       endif
     endif
  endif
  if (present(err)) then

src/setup/setup_grtde.f90

@@ -60,6 +60,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
  use vectorutils,    only:rotatevec
  use gravwaveutils,  only:theta_gw,calc_gravitwaves
  use setup_params,   only:rhozero,npart_total
+ use systemutils,    only:get_command_option
  integer,           intent(in)    :: id
  integer,           intent(inout) :: npart
  integer,           intent(out)   :: npartoftype(:)
@@ -70,7 +71,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
  character(len=20), intent(in)    :: fileprefix
  real,              intent(out)   :: vxyzu(:,:)
  character(len=120) :: filename
- integer :: ierr
+ integer :: ierr,np_default
  logical :: iexist,write_profile,use_var_comp
  real    :: rtidal,rp,semia,period,hacc1,hacc2
  real    :: vxyzstar(3),xyzstar(3)
@@ -100,7 +101,8 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
  call set_units(mass=mhole*solarm,c=1.d0,G=1.d0) !--Set central mass to M=1 in code units
  star%mstar      = 1.*solarm/umass
  star%rstar      = 1.*solarr/udist
- star%np         = 1e6
+ np_default      = 1e6
+ star%np         = int(get_command_option('np',default=np_default)) ! can set default value with --np=1e5 flag (mainly for testsuite)
  star%iprofile   = 2
  beta            = 5.
  ecc             = 0.8
@@ -109,7 +111,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
  theta           = 0.
  write_profile   = .false.
  use_var_comp    = .false.
- relax           = .false.
+ relax           = .true.
 !
 !-- Read runtime parameters from setup file
 !