Merge pull request #413 from danieljprice/set_star

AUTHORS

@@ -20,8 +20,8 @@ Rebecca Nealon <rebecca.nealon@warwick.ac.uk>
 Ward Homan <ward.homan@kuleuven.be>
 Christophe Pinte <christophe.pinte@univ-grenoble-alpes.fr>
 Elisabeth Borchert <elisabeth.borchert@monash.edu>
-Fangyi (Fitz) Hu <fhuu0005@student.monash.edu>
 Megha Sharma <msha0023@student.monash.edu>
+Fangyi (Fitz) Hu <fhuu0005@student.monash.edu>
 Terrence Tricco <ttricco@cita.utoronto.ca>
 Mats Esseldeurs <matsesseldeurs@yahoo.com>
 Simone Ceppi <simone.ceppi@unimi.it>

build/Makefile

@@ -204,9 +204,6 @@ endif
 
 ifeq ($(GR), yes)
     FPPFLAGS += -DGR
-ifeq ($(METRIC), kerr)
-    FPPFLAGS += -DKERR
-endif
 ifeq ($(ISENTROPIC), yes)
     FPPFLAGS += -DISENTROPIC
     FPPFLAGS += -DLIGHTCURVE

build/Makefile_setups

@@ -264,7 +264,6 @@ endif
 
 ifeq ($(SETUP), prtest)
 #   simple test of prdrag
-    FPPFLAGS=
     SETUPFILE= setup_prtest.f90
     KNOWN_SETUP=yes
 endif
@@ -318,15 +317,13 @@ endif
 
 ifeq ($(SETUP), torus)
 #   MRI torus
-    FPPFLAGS=
     SETUPFILE= setup_torus.f90
     ANALYSIS= analysis_torus.f90
     KNOWN_SETUP=yes
 endif
 
 ifeq ($(SETUP), galcen)
 #   galactic centre
-    FPPFLAGS=
     SETUPFILE= setup_galcen_stars.f90
     SRCINJECT= inject_galcen_winds.f90
     KNOWN_SETUP=yes
@@ -689,7 +686,6 @@ endif
 
 ifeq ($(SETUP), cluster)
 #   star cluster formation
-    FPPFLAGS=
     SETUPFILE= velfield_fromcubes.f90 setup_cluster.f90
     MODFILE= moddump_default.f90
     ANALYSIS= phantom_pdfs.o analysis_MWpdf.f90 #analysis_sinkmass.f90
@@ -703,7 +699,6 @@ endif
 
 ifeq ($(SETUP), binary)
 #   binary stars
-    FPPFLAGS=
     #SRCINJECT= utils_binary.f90 set_binary.f90 inject_rochelobe.f90
     SETUPFILE= setup_binary.f90
     #SETUPFILE= setup_chinchen.f90
@@ -713,7 +708,6 @@ endif
 
 ifeq ($(SETUP), hierarchical)
 #   hierarchical system setup
-    FPPFLAGS= -DCONST_AV
     #SRCINJECT= set_hierarchical.f90 set_binary.f90 inject_rochelobe.f90 utils_binary.f90
     SETUPFILE= setup_hierarchical.f90
     KNOWN_SETUP=yes
@@ -722,14 +716,12 @@ endif
 
 ifeq ($(SETUP), common)
 #   binary setup
-    FPPFLAGS=
     SETUPFILE= setup_common.f90
     KNOWN_SETUP=yes
 endif
 
 ifeq ($(SETUP), star)
 #   import stellar model from 1D stellar evolution code
-    FPPFLAGS=
     SETUPFILE= setup_star.f90
     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
@@ -742,7 +734,6 @@ ifeq ($(SETUP), grstar)
 #   star in GR using Minkowski metric
     GR=yes
     METRIC=minkowski
-    FPPFLAGS=
     IND_TIMESTEPS=yes
     SETUPFILE= setup_star.f90
     MODFILE= moddump_tidal.f90
@@ -1003,7 +994,6 @@ endif
 ifeq ($(SETUP), testgr)
 #   unit tests of general relativistic code
     GR=yes
-    FPPFLAGS= -DGR
     KNOWN_SETUP=yes
     METRIC=kerr
     SETUPFILE= setup_grdisc.f90

src/main/checksetup.F90

@@ -804,10 +804,6 @@ subroutine check_gr(npart,nerror,xyzh,vxyzu)
        call pack_metric(xyzh(1:3,i),metrici)
        call unpack_metric(metrici,gcov=gcov)
        call get_u0(gcov,vxyzu(1:3,i),U0,ierr)
-       if (ierr /= 0) then
-          print*,vxyzu(1:3,i),gcov,U0
-          read*
-       endif
        if (ierr/=0) nbad = nbad + 1
     endif
  enddo

src/main/externalforces_gr.F90

@@ -16,11 +16,11 @@ module externalforces
 !   - accradius1      : *soft accretion radius of black hole*
 !   - accradius1_hard : *hard accretion radius of black hole*
 !
-! :Dependencies: dump_utils, infile_utils, io, metric_tools, part, units
+! :Dependencies: dump_utils, infile_utils, io, metric, metric_tools, part,
+!   units
 !
+ use metric, only:mass1
  implicit none
- character(len=80), parameter, public :: &  ! module version
-    modid="$Id$"
 
  private
  public :: externalforce,externalforce_vdependent
@@ -36,7 +36,7 @@ module externalforces
  !
  integer, parameter, public :: iext_gr = 1
 
- real, public :: mass1 = 1.0
+ public :: mass1  ! exported from metric module
  real, public :: accradius1 = 0.
  real, public :: accradius1_hard = 0.
 

src/main/metric_schwarzschild.f90

@@ -22,6 +22,7 @@ module metric
  integer,          parameter :: imetric     = 2
 
  real, public :: mass1 = 1.       ! mass of central object
+ real, public :: a     = 0.       ! spin of central object
 
 contains
 

src/main/options.f90

@@ -127,7 +127,7 @@ subroutine set_default_options
  ! artificial thermal conductivity
  alphau = 1.
  if (gr) alphau = 0.1
- ireconav = 1
+ ireconav = -1
 
  ! artificial resistivity (MHD only)
  alphaB            = 1.0

src/setup/relax_star.f90

@@ -18,10 +18,10 @@ module relaxstar
 !   - tol_dens : *% error in density to stop relaxation*
 !   - tol_ekin : *tolerance on ekin/epot to stop relaxation*
 !
-! :Dependencies: checksetup, damping, deriv, dim, energies, eos, fileutils,
-!   infile_utils, initial, io, io_summary, memory, options, part, physcon,
-!   ptmass, readwrite_dumps, setstar_utils, sortutils, step_lf_global,
-!   table_utils, units
+! :Dependencies: checksetup, damping, deriv, dim, energies, eos,
+!   externalforces, fileutils, infile_utils, initial, io, io_summary,
+!   memory, options, part, physcon, ptmass, readwrite_dumps, setstar_utils,
+!   sortutils, step_lf_global, table_utils, units
 !
  implicit none
  public :: relax_star,write_options_relax,read_options_relax
@@ -30,7 +30,7 @@ module relaxstar
  real,    private :: tol_dens = 1.   ! allow 1% RMS error in density
  integer, private :: maxits = 1000
 
- real,    private :: gammaprev,hfactprev
+ real,    private :: gammaprev,hfactprev,mass1prev
  integer, private :: ieos_prev
 
  integer, public :: ierr_setup_errors = 1, &
@@ -113,7 +113,7 @@ subroutine relax_star(nt,rho,pr,r,npart,xyzh,use_var_comp,Xfrac,Yfrac,mu,ierr,np
  !
  ! check particle setup is sensible
  !
- call check_setup(nwarn,nerr,restart=.true.) ! restart=T allows accreted/masked particles
+ call check_setup(nerr,nwarn,restart=.true.) ! restart=T allows accreted/masked particles
  if (nerr > 0) then
     call error('relax_star','cannot relax star because particle setup contains errors')
     call restore_original_options(i1,npart)
@@ -372,15 +372,17 @@ end subroutine reset_u_and_get_errors
 !----------------------------------------------------------------
 subroutine set_options_for_relaxation(tdyn)
  use eos,  only:ieos,gamma
- use part, only:hfact,maxvxyzu
+ use part, only:hfact,maxvxyzu,gr
  use damping, only:damp,tdyn_s
  use options, only:idamp
- use units,   only:utime
+ use units,          only:utime
+ use externalforces, only:mass1
  real, intent(in) :: tdyn
 
  gammaprev = gamma
  hfactprev = hfact
  ieos_prev = ieos
+ mass1prev = mass1
  !
  ! turn on settings appropriate to relaxation
  !
@@ -392,6 +394,7 @@ subroutine set_options_for_relaxation(tdyn)
     idamp = 1
     damp = 0.05
  endif
+ if (gr) mass1 = 0. ! use Minkowski metric during relaxation
 
 end subroutine set_options_for_relaxation
 
@@ -422,7 +425,8 @@ subroutine restore_original_options(i1,npart)
  use eos,     only:ieos,gamma
  use damping, only:damp
  use options, only:idamp
- use part,    only:hfact,vxyzu
+ use part,    only:hfact,vxyzu,gr
+ use externalforces, only:mass1
  integer, intent(in) :: i1,npart
 
  gamma = gammaprev
@@ -431,6 +435,7 @@ subroutine restore_original_options(i1,npart)
  idamp = 0
  damp = 0.
  vxyzu(1:3,i1+1:npart) = 0.
+ if (gr) mass1 = mass1prev
 
 end subroutine restore_original_options
 

src/setup/set_binary.f90

@@ -91,7 +91,7 @@ subroutine set_binary(m1,m2,semimajoraxis,eccentricity, &
 
  if (do_verbose) then
     print "(/,2x,a)",'---------- binary parameters ----------- '
-    print "(8(2x,a,g14.6,/),2x,a,g14.6)", &
+    print "(8(2x,a,1pg14.6,/),2x,a,1pg14.6)", &
         'primary mass     :',m1, &
         'secondary mass   :',m2, &
         'mass ratio m2/m1 :',m2/m1, &
@@ -167,7 +167,7 @@ subroutine set_binary(m1,m2,semimajoraxis,eccentricity, &
     E_dot = sqrt((m1 + m2)/(a**3))/(1.-eccentricity*cos(E))
 
     if (do_verbose) then
-       print "(4(2x,a,g12.4,/),2x,a,g12.4)", &
+       print "(4(2x,a,1pg14.6,/),2x,a,1pg14.6)", &
              'Eccentric anomaly:',E, &
              'E_dot            :',E_dot, &
              'inclination     (i, deg):',incl, &
@@ -202,7 +202,7 @@ subroutine set_binary(m1,m2,semimajoraxis,eccentricity, &
 
  ! print info about positions and velocities
  if (do_verbose) then
-    print "(7(2x,a,g12.4,/),2x,a,g12.4)", &
+    print "(7(2x,a,1pg14.6,/),2x,a,1pg14.6)", &
         'angular momentum :',angmbin, &
         'mean ang. speed  :',omega0, &
         'Omega_0 (prim)   :',v1(2)/x1(1), &
@@ -218,7 +218,7 @@ subroutine set_binary(m1,m2,semimajoraxis,eccentricity, &
     omega_corotate = omega0
     v1(2) = v1(2) - omega0*x1(1)
     v2(2) = v2(2) - omega0*x2(1)
-    if (do_verbose) print "(2(2x,a,g12.4,/))", &
+    if (do_verbose) print "(2(2x,a,1pg14.6,/))", &
      'Omega_0 (primary)     :',v1(2)/x1(1), &
      'Omega_0 (secondary)   :',v2(2)/x2(1)
  endif