Skip to content
Permalink
Browse files

Merge pull request #15 from evanoconnor/smallfixes

Performance and analysis improvements
  • Loading branch information...
evanoconnor committed Feb 15, 2019
2 parents 1677b0b + 19b6b83 commit 59f6eee6df036e45f807e4fdf3ebde8279545535
Showing with 24 additions and 9 deletions.
  1. +1 −0 src/GR1D_module.F90
  2. +12 −0 src/M1/M1_conservativeupdate.F90
  3. +8 −8 src/M1/M1_explicitterms.F90
  4. +1 −0 src/initialize_vars.F90
  5. +2 −1 src/output.F90
@@ -173,6 +173,7 @@ module GR1D_module
real*8 :: total_energy_radiated
real*8 :: total_energy_absorped
real*8 :: total_net_heating
real*8 :: total_net_deintde
real*8 :: total_mass_gain

!M1 controls
@@ -18,9 +18,11 @@ subroutine M1_conservativeupdate(dts)
integer keyerr,keytemp
real*8 eosdummy(14)
logical :: passfluxtest
real*8 :: oldeps_forheat(n1)

nogain = .true.
total_net_heating = 0.0d0
total_net_deintdt = 0.0d0
total_mass_gain = 0.0d0
igain(1) = ghosts1+1
gain_radius = 0.0d0
@@ -93,6 +95,8 @@ subroutine M1_conservativeupdate(dts)
endif

!reconstruct primatives
!but first, retain old eps to determined heating rate
oldeps_forheat = eps
call con2prim

! eos update, eps fixed, find temp,entropy,cs2 etc.
@@ -117,6 +121,14 @@ subroutine M1_conservativeupdate(dts)
write(6,*) "timestep number: ",nt
write(6,"(i4,1P10E15.6)") k,x1(k),rho(k)/rho_gf,temp(k),eps(k)/eps_gf,ye(k)
endif

if (eps(k).gt.oldeps_forheat(k)) then
if ((rho(k)/rho_gf.lt.3.0d10).and.(entropy(k).gt.6.0d0)) then
!gain region
total_net_deintdt = total_net_deintdt + volume(k)*rho(k)*(eps(k)-oldeps_forheat(k))/energy_gf/(dts/time_gf)
endif
endif

enddo

!GR gravity updates that rely on primitive variables
@@ -211,11 +211,11 @@ subroutine M1_explicitterms(dts,implicit_factor)
endif

!actual speed
l_min(1) = (3.0d0*M1chi_space_minus(k+1)-1.0d0)*0.5d0*l_min_thin(1) + &
(1.0d0 - M1chi_space_minus(k+1))*1.5d0*l_min_thick(1)
l_min(1) = (3.0d0*M1chi_space(k+1)-1.0d0)*0.5d0*l_min_thin(1) + &
(1.0d0 - M1chi_space(k+1))*1.5d0*l_min_thick(1)

l_max(1) = (3.0d0*M1chi_space_minus(k+1)-1.0d0)*0.5d0*l_max_thin(1) + &
(1.0d0 - M1chi_space_minus(k+1))*1.5d0*l_max_thick(1)
l_max(1) = (3.0d0*M1chi_space(k+1)-1.0d0)*0.5d0*l_max_thin(1) + &
(1.0d0 - M1chi_space(k+1))*1.5d0*l_max_thick(1)

!plus interface (k zone)
!thin limit:
@@ -274,10 +274,10 @@ subroutine M1_explicitterms(dts,implicit_factor)
endif

!actualy speed
l_min(2) = (3.0d0*M1chi_space_plus(k)-1.0d0)*0.5d0*l_min_thin(2) + &
(1.0d0 - M1chi_space_plus(k))*1.5d0*l_min_thick(2)
l_max(2) = (3.0d0*M1chi_space_plus(k)-1.0d0)*0.5d0*l_max_thin(2) + &
(1.0d0 - M1chi_space_plus(k))*1.5d0*l_max_thick(2)
l_min(2) = (3.0d0*M1chi_space(k)-1.0d0)*0.5d0*l_min_thin(2) + &
(1.0d0 - M1chi_space(k))*1.5d0*l_min_thick(2)
l_max(2) = (3.0d0*M1chi_space(k)-1.0d0)*0.5d0*l_max_thin(2) + &
(1.0d0 - M1chi_space(k))*1.5d0*l_max_thick(2)


!check for NaNs
@@ -61,6 +61,7 @@ subroutine initialize_vars
total_energy_radiated = 0.0d0
total_energy_absorped = 0.0d0
total_net_heating = 0.0d0
total_net_deintdt = 0.0d0

shock_radius = 0.0d0
ishock(1) = ghosts1+1
@@ -521,11 +521,12 @@ subroutine output_all(modeflag)


scalars(1:nscalars0) = 0.0d0
nscalars = 4
nscalars = 5
scalars(1) = total_net_heating/(luminosity(1)+luminosity(2))
scalars(2) = luminosity(1)+luminosity(2)
scalars(3) = total_net_heating
scalars(4) = total_mass_gain
scalars(5) = total_net_deintdt
filename = trim(adjustl(outdir))//"/M1_net_heating.dat"
call output_many_scalars(scalars,nscalars0,nscalars,filename)

0 comments on commit 59f6eee

Please sign in to comment.
You can’t perform that action at this time.