Merge pull request #24 from themikelau/master

changes to star setup and relaxation

src/main/eos.F90

@@ -1116,31 +1116,25 @@ end subroutine calc_rec_ene
 !  pressure and density
 !+
 !----------------------------------------------------------------
-subroutine calc_temp_and_ene(ieos,mu,XX,YY,gamma,rho,pres,ene,temp,ierr)
+subroutine calc_temp_and_ene(rho,pres,ene,temp,ierr,guesseint)
  use physcon,          only:kb_on_mh
  use eos_idealplusrad, only:get_idealgasplusrad_tempfrompres,get_idealplusrad_enfromtemp
- real, intent(in)      :: rho,pres,mu,gamma,XX,YY
- real, intent(inout)   :: ene,temp
- integer, intent(in)   :: ieos
- integer, intent(out)  :: ierr
- real                  :: e_rec
-
+ use eos_mesa,         only:get_eos_eT_from_rhop_mesa
+ real, intent(in)           :: rho,pres
+ real, intent(inout)        :: ene,temp
+ real, intent(in), optional :: guesseint
+ integer, intent(out)       :: ierr
+ 
  ierr = 0
  select case(ieos)
  case(2) ! Adiabatic/polytropic EoS
-    temp = pres / (rho * kb_on_mh) * mu
+    temp = pres / (rho * kb_on_mh) * gmw
     ene = pres / ( (gamma-1.) * rho)
  case(12) ! Ideal plus rad. EoS
-    call get_idealgasplusrad_tempfrompres(pres,rho,mu,temp)
-    call get_idealplusrad_enfromtemp(rho,temp,mu,ene)
+    call get_idealgasplusrad_tempfrompres(pres,rho,gmw,temp)
+    call get_idealplusrad_enfromtemp(rho,temp,gmw,ene)
  case(10) ! MESA-like EoS
-    ! Approximate the temperature as that from ideal gas plus radiation
-    call get_idealgasplusrad_tempfrompres(pres,rho,mu,temp)
-
-    ! Calculate internal energy from gas and radiation, then add recombination energy
-    call get_idealplusrad_enfromtemp(rho,temp,mu,ene)
-    call calc_rec_ene(XX,YY,e_rec)
-    ene = ene + e_rec
+    call get_eos_eT_from_rhop_mesa(rho,pres,ene,temp,guesseint)
  case default
     ierr = 1
  end select

src/main/eos_mesa.f90

@@ -88,7 +88,6 @@ subroutine get_eos_pressure_gamma1_mesa(den,eint,pres,gam1,ierr)
  integer, intent(out) :: ierr
 
  call getvalue_mesa(den,eint,2,pres,ierr)
- !call getvalue_mesa(den,eint,3,pres,ierr)
  call getvalue_mesa(den,eint,11,gam1)
 
 end subroutine get_eos_pressure_gamma1_mesa
@@ -109,7 +108,7 @@ end subroutine get_eos_kappa_mesa
 
 !----------------------------------------------------------------
 !+
-!  subroutine returns pressure/temperature as
+!  subroutine returns pressure and temperature as
 !  a function of density/internal energy
 !+
 !----------------------------------------------------------------
@@ -122,6 +121,74 @@ subroutine get_eos_pressure_temp_mesa(den,eint,pres,temp)
 
 end subroutine get_eos_pressure_temp_mesa
 
+!----------------------------------------------------------------
+!+
+!  subroutine returns internal energy and temperature from
+!  density and internal energy using bisection method
+!+
+!----------------------------------------------------------------
+subroutine get_eos_eT_from_rhop_mesa(rho,pres,eint,temp,guesseint)
+ real, intent(in)           :: rho,pres
+ real, intent(out)          :: eint,temp
+ real, intent(in), optional :: guesseint
+ real                       :: err,eintguess,eint1,eint2,&
+                               eint3,pres1,pres2,pres3,left,right,mid
+ real, parameter            :: tolerance = 1d-15
+ integer                    :: ierr
+
+ if (present(guesseint)) then
+    eintguess = guesseint
+    eint1 = 1.005 * eintguess  ! Tight lower bound
+    eint2 = 0.995 * eintguess  ! Tight upper bound
+ else
+    eintguess = 1.5*pres/rho
+    eint1 = 10. * eintguess  ! Guess lower bound
+    eint2 = 0.1 * eintguess  ! Guess upper bound
+ endif
+
+ call getvalue_mesa(rho,eint1,2,pres1,ierr)
+ call getvalue_mesa(rho,eint2,2,pres2,ierr)
+ left  = pres - pres1
+ right = pres - pres2
+
+ ! If lower and upper bounds do not contain roots, extend them until they do
+ do while (left*right > 0.)
+    eint1 = 0.99 * eint1
+    eint2 = 1.01 * eint2
+    call getvalue_mesa(rho,eint1,2,pres1,ierr)
+    call getvalue_mesa(rho,eint2,2,pres2,ierr)
+    left  = pres - pres1
+    right = pres - pres2
+ enddo
+
+ ! Start bisecting
+ err = huge(1.)
+ do while (abs(err) > tolerance)
+    call getvalue_mesa(rho,eint1,2,pres1,ierr)
+    call getvalue_mesa(rho,eint2,2,pres2,ierr)
+    left  = pres - pres1
+    right = pres - pres2
+    eint3 = 0.5*(eint2+eint1)
+    call getvalue_mesa(rho,eint3,2,pres3,ierr)
+    mid = pres - pres3
+
+    if (left*mid < 0.) then
+       eint2 = eint3
+    elseif (right*mid < 0.) then
+       eint1 = eint3
+    elseif (mid == 0.) then
+       eint = eint3
+       exit 
+    endif
+
+    eint = eint3
+    err = (eint2 - eint1)/eint1
+ enddo
+
+ call getvalue_mesa(rho,eint,4,temp,ierr)
+
+end subroutine get_eos_eT_from_rhop_mesa 
+
 !----------------------------------------------------------------
 !+
 !  subroutine returns various quantities as

src/main/eos_mesa_microphysics.F90

@@ -429,6 +429,7 @@ end subroutine read_eos_mesa
 ! 1. logRho        2. logP          3. logPgas       4. logT
 ! 5. dlnP/dlnrho|e 6. dlnP/dlne|rho 7. dlnT/dlnrho|e 8. dlnT/dlne|rho
 ! 9. logS         10. dlnT/dlnP|S  11. Gamma1       12. gamma
+! Note: ivout=1,2,3,4 returns the unlogged quantity
 subroutine getvalue_mesa(rho,eint,ivout,vout,ierr)
  real, intent(in) :: rho, eint
  real, intent(out) :: vout

src/main/extern_densprofile.f90

@@ -32,7 +32,7 @@ module extern_densprofile
 
  ! *** Add option to .in file to specify density profile / mass enclosed filename? ***
  character(len=*), public, parameter     :: rhotabfile = 'density-profile.tab'
- integer, public, parameter              :: nrhotab = 1001  ! maximum allowed size of r rho tabulated arrays
+ integer, public, parameter              :: nrhotab = 5000  ! maximum allowed size of r rho tabulated arrays
 
  public :: densityprofile_force, load_extern_densityprofile, read_rhotab, write_rhotab, calc_menc
  public :: read_rhotab_wrapper

src/main/readwrite_infile.F90

@@ -210,7 +210,7 @@ subroutine write_infile(infile,logfile,evfile,dumpfile,iwritein,iprint)
  ! thermodynamics
  !
  call write_options_eos(iwritein)
- if (maxvxyzu >= 4 .and. (ieos==2 .or. ieos==10 .or. ieos==15) ) then
+ if (maxvxyzu >= 4 .and. (ieos==2 .or. ieos==10 .or. ieos==15 .or. ieos==12 .or. ieos==16) ) then
     call write_inopt(ipdv_heating,'ipdv_heating','heating from PdV work (0=off, 1=on)',iwritein)
     call write_inopt(ishock_heating,'ishock_heating','shock heating (0=off, 1=on)',iwritein)
     if (mhd) then

src/setup/density_profiles.f90

@@ -512,15 +512,37 @@ end subroutine read_mesa
 !  used in star setup to write softened stellar profile.
 !+
 !----------------------------------------------------------------
-subroutine write_softened_profile(outputpath, m, pres, temp, r, rho, ene)
- real, allocatable               :: m(:),rho(:),pres(:),r(:),ene(:),temp(:)
+subroutine write_softened_profile(outputpath, m, pres, temp, r, rho, ene, Xfrac, Yfrac, csound)
+ real, intent(in)                :: m(:),rho(:),pres(:),r(:),ene(:),temp(:)
+ real, intent(in), optional      :: Xfrac(:),Yfrac(:),csound(:)
  character(len=120), intent(in)  :: outputpath
  integer                         :: i
+
  open(1, file = outputpath, status = 'new')
- write(1,'(a)') '[    Mass   ]  [  Pressure ]  [Temperature]  [   Radius  ]  [  Density  ]  [   E_int   ]'
- write(1,42) (m(i), pres(i), temp(i), r(i), rho(i), ene(i), i = 1, size(r))
-42 format (es13.7, 2x, es13.7, 2x, es13.7, 2x, es13.7, 2x, es13.7, 2x, es13.7)
+
+ if (present(Xfrac) .and. present(Yfrac)) then
+    if (present(csound)) then
+       write(1,'(a)') '[    Mass   ]  [  Pressure ]  [Temperature]  [   Radius  ]  &
+                      &[  Density  ]  [   E_int   ]  [   Xfrac   ]  [   Yfrac   ]  [Sound speed]'
+       write(1,101) (m(i),pres(i),temp(i),r(i),rho(i),ene(i),Xfrac(i),Yfrac(i),csound(i),i=1,size(r))
+       101 format (es13.7,2x,es13.7,2x,es13.7,2x,es13.7,2x,es13.7,2x,es13.7,2x,es13.7,&
+       2x,es13.7,2x,es13.7)
+    else
+       write(1,'(a)') '[    Mass   ]  [  Pressure ]  [Temperature]  [   Radius  ]  &
+                      &[  Density  ]  [   E_int   ]  [   Xfrac   ]  [   Yfrac   ]'
+       write(1,102) (m(i),pres(i),temp(i),r(i),rho(i),ene(i),Xfrac(i),Yfrac(i),i=1,size(r))
+       102 format (es13.7,2x,es13.7,2x,es13.7,2x,es13.7,2x,es13.7,2x,es13.7,2x,es13.7,&
+       2x,es13.7)
+    endif
+ else
+    write(1,'(a)') '[    Mass   ]  [  Pressure ]  [Temperature]  [   Radius  ]  &
+                   &[  Density  ]  [   E_int   ]'
+    write(1,103) (m(i),pres(i),temp(i),r(i),rho(i),ene(i),i=1,size(r))
+    103 format (es13.7,2x,es13.7,2x,es13.7,2x,es13.7,2x,es13.7,2x,es13.7)
+ endif
+
  close(1, status = 'keep')
+
 end subroutine write_softened_profile
 
 !-----------------------------------------------------------------------

src/setup/relax_star.f90

@@ -230,8 +230,8 @@ subroutine shift_particles(npart,xyzh,vxyzu,dtmin)
  integer, intent(in) :: npart
  real, intent(inout) :: xyzh(:,:), vxyzu(:,:)
  real, intent(out)   :: dtmin
- real :: dx(3),dti,phi,rhoi,cs
- integer :: i
+ real :: dx(3),dti,phi,rhoi,cs,hi
+ integer :: i,nlargeshift
 !
 ! get forces on particles
 !
@@ -240,25 +240,33 @@ subroutine shift_particles(npart,xyzh,vxyzu,dtmin)
 ! shift particles asynchronously
 !
  dtmin = huge(dtmin)
+ nlargeshift = 0
  !$omp parallel do schedule(guided) default(none) &
  !$omp shared(npart,xyzh,vxyzu,fxyzu,fext,xyzmh_ptmass,nptmass,massoftype,ieos) &
- !$omp private(i,dx,dti,phi,cs,rhoi) &
- !$omp reduction(min:dtmin)
+ !$omp private(i,dx,dti,phi,cs,rhoi,hi) &
+ !$omp reduction(min:dtmin) &
+ !$omp reduction(+:nlargeshift)
  do i=1,npart
     fext(1:3,i) = 0.
     if (nptmass > 0) then
        call get_accel_sink_gas(nptmass,xyzh(1,i),xyzh(2,i),xyzh(3,i),xyzh(4,i),&
                                xyzmh_ptmass,fext(1,i),fext(2,i),fext(3,i),phi)
     endif
-    rhoi = rhoh(xyzh(4,i),massoftype(igas))
+    hi = xyzh(4,i)
+    rhoi = rhoh(hi,massoftype(igas))
     cs = get_spsound(ieos,xyzh(:,i),rhoi,vxyzu(:,i))
-    dti = 0.3*xyzh(4,i)/cs   ! h/cs
+    dti = 0.3*hi/cs   ! h/cs
     dx  = 0.5*dti**2*(fxyzu(1:3,i) + fext(1:3,i))
+    if (dot_product(dx,dx) > hi**2) then
+       dx = dx / sqrt(dot_product(dx,dx)) * hi  ! Avoid large shift in particle position
+       nlargeshift = nlargeshift + 1
+    endif
     xyzh(1:3,i) = xyzh(1:3,i) + dx(:)
     vxyzu(1:3,i) = dx(:)/dti ! fake velocities, so we can measure kinetic energy
     dtmin = min(dtmin,dti)   ! used to print a "time" in the output (but it is fake)
  enddo
  !$omp end parallel do
+ if (nlargeshift > 0) print*,'Warning: Restricted dx for ', nlargeshift, 'particles'
 
 end subroutine shift_particles
 

src/setup/set_softened_core.f90

@@ -26,9 +26,7 @@
 !--------------------------------------------------------------------------
 module setsoftenedcore
  use physcon,          only:pi,gg,solarm,solarr,kb_on_mh
- use eos_idealplusrad, only:get_idealgasplusrad_tempfrompres,&
-                            get_idealplusrad_enfromtemp
- use table_utils,      only:interpolator,diff,flip_array
+ use table_utils,      only:interpolator,diff
 
  implicit none
  real    :: hsoft,msoft,mcore
@@ -48,25 +46,19 @@ module setsoftenedcore
 
 !-----------------------------------------------------------------------
 !+
-!  Main subroutine that calls the subroutines to calculate the cubic
-!  core profile
+!  Main subroutine that calculates the cubic core profile
 !+
 !-----------------------------------------------------------------------
-subroutine set_softened_core(ieos,gamma,constX,constY,mu,mcore,hsoft,hphi,&
-                             rho,r,pres,m,ene,temp,ierr,Xfrac,Yfrac)
+subroutine set_softened_core(mcore,hsoft,hphi,rho,r,pres,m,ene,temp,ierr)
+ use eos,         only:calc_temp_and_ene
+ use table_utils, only:flip_array
  real, intent(inout)        :: r(:),rho(:),m(:),pres(:),ene(:),temp(:)
  real, allocatable          :: phi(:)
- real, intent(in), optional :: Xfrac(:),Yfrac(:)
- real, intent(in)           :: mu,mcore,hsoft,hphi,gamma,constX,constY
- integer, intent(in)        :: ieos
+ real, intent(in)           :: mcore,hsoft,hphi
  integer, intent(out)       :: ierr
- real                       :: mc,h,hphi_cm
+ real                       :: mc,h,hphi_cm,eneguess
  logical                    :: isort_decreasing,iexclude_core_mass
-
- ! Xfrac, Yfrac:    H and He mass fractions from the MESA profile
- ! constX, constY:  H and He mass fractions at r = R/2
- ! mu:              Mean molecular weight at r = R/2
- ! gamma:           Adiabatic exponent for adiabatic EoS
+ integer                    :: i
 
  ! Output data to be sorted from stellar surface to interior?
  isort_decreasing = .true.     ! Needs to be true if to be read by Phantom
@@ -84,11 +76,13 @@ subroutine set_softened_core(ieos,gamma,constX,constY,mu,mcore,hsoft,hphi,&
  call calc_phi(r, mc, m-mc, hphi_cm, phi)
  call calc_pres(r, rho, phi, pres)
 
- if (present(Xfrac) .and. present(Yfrac)) then ! This case is temporarily forbidden until variable composition is implemented in Phantom
-    call calc_softened_temp_and_ene(ieos,hidx,mu,constX,constY,gamma,rho,pres,ene,temp,ierr,Xfrac,Yfrac)
- else
-    call calc_softened_temp_and_ene(ieos,hidx,mu,constX,constY,gamma,rho,pres,ene,temp,ierr)
- endif
+ call calc_temp_and_ene(rho(1),pres(1),ene(1),temp(1),ierr)
+ do i = 2,size(rho)-1
+    eneguess = ene(i-1)
+    call calc_temp_and_ene(rho(i),pres(i),ene(i),temp(i),ierr,eneguess)
+ enddo
+ ene(size(rho))  = 0. ! Zero surface internal energy
+ temp(size(rho)) = 0. ! Zero surface temperature
 
  ! Reverse arrays so that data is sorted from stellar surface to stellar centre.
  if (isort_decreasing) then
@@ -305,46 +299,4 @@ subroutine calc_pres(r, rho, phi, pres)
  enddo
 end subroutine calc_pres
 
-!----------------------------------------------------------------
-!+
-!  Calculates temperature and specific internal energy for the
-!  softened star from pressure and density
-!+
-!----------------------------------------------------------------
-subroutine calc_softened_temp_and_ene(ieos,hidx,mu,constX,constY,gamma,rho,pres,ene,temp,ierr,Xfrac,Yfrac)
- use eos,                only:calc_temp_and_ene
- real, intent(in)           :: rho(:),pres(:),mu,gamma,constX,constY
- real, intent(in), optional :: Xfrac(:),Yfrac(:)
- real, intent(inout)        :: ene(:),temp(:)
- integer, intent(in)        :: ieos,hidx
- integer, intent(out)       :: ierr
- real                       :: muprofile(size(rho))
- integer                    :: endidx,i
-
- ierr = 0
- if (ieos == 10) then
-    endidx = hidx ! For r > h, we just use the original MESA profile
- else
-    endidx = size(rho) ! Recalculate u and T for r > h too
- endif
-
- if (ieos == 12) temp(size(rho)) = 0. ! Zero surface temperature
-
- if ( (.not. present(Xfrac)) .and. (.not. present(Yfrac)) ) then
-    ! Calculate T and u for the entire star using fixed mu, X, Y
-    do i = 1,endidx
-       call calc_temp_and_ene(ieos,mu,constX,constY,gamma,rho(i),pres(i),ene(i),temp(i),ierr)
-    enddo
- elseif (present(Xfrac) .and. present(Yfrac)) then ! This case is temporarily forbidden until particles with variable composition is implemented in Phantom
-    ! Calculate mean molecular weight profile from X, Y profiles
-    muprofile = 4. / (6.*Xfrac + Yfrac + 2.)
-    do i = 1,endidx ! For r > h, we just use the original MESA profile
-       call calc_temp_and_ene(ieos,muprofile(i),Xfrac(i),Yfrac(i),gamma,rho(i),pres(i),ene(i),temp(i),ierr)
-    enddo
- else
-    ierr = 1
- endif
-
-end subroutine calc_softened_temp_and_ene
-
 end module setsoftenedcore

src/setup/setup_star.f90

@@ -310,9 +310,9 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
        hsoft = 0.5*hdens ! This is set by default so that the pressure, energy, and temperature
        ! are same as the original profile for r > hsoft
 
-       call set_softened_core(ieos,gamma,X_in,Y_in,gmw,mcore,hdens,hsoft,rho0,r0,pres0,m0,ene0,temp0,ierr)
+       call init_eos(ieos,ierr)
+       call set_softened_core(mcore,hdens,hsoft,rho0,r0,pres0,m0,ene0,temp0,ierr)
        if (ierr==1) call fatal('setup','EoS not one of: adiabatic, ideal gas plus radiation, MESA in set_softened_core')
-       !if (ierr==2) call fatal('setup','Xfrac and Yfrac not provided to set_softened_core for ieos=10 (MESA EoS)')
        call set_stellar_core(nptmass,xyzmh_ptmass,vxyz_ptmass,mcore,hsoft,ihsoft)
        call write_softened_profile(outputfilename,m0,pres0,temp0,r0,rho0,ene0)
        densityfile = outputfilename ! Have the read_mesa_file subroutine read the softened profile instead
@@ -372,7 +372,6 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
  !
  ! relax the density profile to achieve nice hydrostatic equilibrium
  !
- call init_eos(ieos,ierr)
  if (relax_star_in_setup) then
     if (nstar==npart) then
        call relax_star(npts,den,pres,r,npart,xyzh)