Skip to content

Commit

Permalink
final updates to GREP in GR1D, small other changes including better g…
Browse files Browse the repository at this point in the history 
…ain region tracking and reporting
  • Loading branch information
Evan O'Connor committed Nov 5, 2015
1 parent e4935dc commit dbffc85
Show file tree
Hide file tree
Showing 8 changed files with 153 additions and 61 deletions.
2 changes: 2 additions & 0 deletions src/GR1D_module.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,6 +170,7 @@ module GR1D_module
real*8 :: total_energy_radiated
real*8 :: total_energy_absorped
real*8 :: total_net_heating
real*8 :: total_mass_gain

!M1 controls
integer :: M1_do_backwardfix
Expand Down Expand Up @@ -238,6 +239,7 @@ module GR1D_module
! ye
real*8,allocatable,save :: ye(:),yem(:),yep(:),ye_prev(:)
real*8,allocatable,save :: dyedt_hydro(:),dyedt_neutrino(:)
real*8,allocatable,save :: depsdt(:),dyedt(:)
real*8,allocatable,save :: ynu(:)
!chemical potential
real*8,allocatable,save :: nuchem(:),elechem(:)
Expand Down
34 changes: 15 additions & 19 deletions src/M1/M1_conservativeupdate.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,9 +17,11 @@ subroutine M1_conservativeupdate(dts)
logical nogain
integer keyerr,keytemp
real*8 eosdummy(14)
logical :: passfluxtest

nogain = .true.
total_net_heating = 0.0d0
total_mass_gain = 0.0d0
igain(1) = ghosts1+1
gain_radius = 0.0d0
maxye = 0.0d0
Expand All @@ -38,26 +40,20 @@ subroutine M1_conservativeupdate(dts)
dtau = 4.0d0*pi*dts*M1_matter_source(k,3) !-\alpha^2S^t * 4*pi*dts
dSr = 4.0d0*pi*dts*M1_matter_source(k,2) ! -\alpha X S^r * 4*pi*dts

!careful when determining gain radius, we require 10 zones of heating
if (dtau.gt.0.0d0.and.nogain) then
gaintracker = gaintracker+1
if (gaintracker.eq.10) then
igain(1) = k-9
gain_radius = x1(igain(1))
nogain = .false.
endif
else if (nogain) then
gaintracker = 0
total_net_heating = 0.0d0
else
!have gain region, don't reset gaintracker
endif
depsdt(k) = M1_matter_source(k,3)/rho(k)*4.0d0*pi/eps_gf*time_gf
dyedt(k) = M1_matter_source(k,4)*4.0d0*pi*X(k)/q(k,1)*(amu_cgs*mass_gf)*time_gf

passfluxtest = sum(abs(q_M1(k,1,:,2))/X(k))/sum(q_M1(k,1,:,1)).gt.0.5d0
passfluxtest = passfluxtest.and.(sum(abs(q_M1(k,2,:,2))/X(k))/sum(q_M1(k,2,:,1)).gt.0.5d0)

passfluxtest = rho(k)/rho_gf.lt.3.0d10

if ((dtau.gt.0.0d0).and.(entropy(k).gt.6.0d0).and.passfluxtest) then
total_net_heating = total_net_heating + &
dtau*X(k)*volume(k)/(energy_gf*dts/time_gf)
total_mass_gain = total_mass_gain + &
volume(k)*rho(k)

if (gaintracker.gt.0.or..not.nogain) then
if (dtau.gt.0.0d0) then
total_net_heating = total_net_heating + &
dtau*volume(k)/(energy_gf*dts/time_gf)
endif
endif

total_energy_absorped = total_energy_absorped + &
Expand Down
101 changes: 77 additions & 24 deletions src/M1/M1_explicitterms.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -161,8 +161,13 @@ subroutine M1_explicitterms(dts,implicit_factor)
l_min_thin(1) = -Xm(k+1)*alpm(k+1)
l_max_thin(1) = Xm(k+1)*alpm(k+1)
else
l_min_thin(1) = -1.0d0
l_max_thin(1) = 1.0d0
if (do_effectivepotential) then
l_min_thin(1) = -alpm(k+1)
l_max_thin(1) = alpm(k+1)
else
l_min_thin(1) = -1.0d0
l_max_thin(1) = 1.0d0
endif
endif

!thick limit:
Expand All @@ -174,8 +179,13 @@ subroutine M1_explicitterms(dts,implicit_factor)
l_min_thick(1) = min((2.0d0*Wm(k+1)**2*p - sqrtdiscrim)/(2.0d0*Wm(k+1)**2+1.0d0),p)
l_max_thick(1) = max((2.0d0*Wm(k+1)**2*p + sqrtdiscrim)/(2.0d0*Wm(k+1)**2+1.0d0),p)
else
p = v1m(k+1)
discrim = (2.0d0*oneWm**2+1.0d0)-2.0d0*oneWm**2*p*p
if (do_effectivepotential) then
p = alpm(k+1)*v1m(k+1)
discrim = alpm(k+1)**2*(2.0d0*oneWm**2+1.0d0)-2.0d0*oneWm**2*p*p
else
p = v1m(k+1)
discrim = (2.0d0*oneWm**2+1.0d0)-2.0d0*oneWm**2*p*p
endif
sqrtdiscrim = sqrt(discrim)
l_min_thick(1) = min((2.0d0*oneWm**2*p - sqrtdiscrim)/(2.0d0*oneWm**2+1.0d0),p)
l_max_thick(1) = max((2.0d0*oneWm**2*p + sqrtdiscrim)/(2.0d0*oneWm**2+1.0d0),p)
Expand All @@ -187,7 +197,11 @@ subroutine M1_explicitterms(dts,implicit_factor)
if (GR) then
discrim = alpm(k+1)**2*Xm(k+1)**2*3.0d0
else
discrim = 3.0d0
if (do_effectivepotential) then
discrim = alpm(k+1)**2*3.0d0
else
discrim = 3.0d0
endif
endif
sqrtdiscrim = sqrt(discrim)
l_min_thick(1) = -sqrtdiscrim/3.0d0
Expand All @@ -209,8 +223,13 @@ subroutine M1_explicitterms(dts,implicit_factor)
l_min_thin(2) = -Xp(k)*alpp(k)
l_max_thin(2) = Xp(k)*alpp(k)
else
l_min_thin(2) = -1.0d0
l_max_thin(2) = 1.0d0
if (do_effectivepotential) then
l_min_thin(2) = -alpp(k)
l_max_thin(2) = alpp(k)
else
l_min_thin(2) = -1.0d0
l_max_thin(2) = 1.0d0
endif
endif


Expand All @@ -223,8 +242,13 @@ subroutine M1_explicitterms(dts,implicit_factor)
l_min_thick(2) = min((2.0d0*Wp(k)**2*p - sqrtdiscrim)/(2.0d0*Wp(k)**2+1.0d0),p)
l_max_thick(2) = max((2.0d0*Wp(k)**2*p + sqrtdiscrim)/(2.0d0*Wp(k)**2+1.0d0),p)
else
p = v1p(k)
discrim = (2.0d0*Wp(k)**2+1.0d0)-2.0d0*Wp(k)**2*p*p
if (do_effectivepotential) then
p = alpp(k)*v1p(k)
discrim = alpp(k)**2*(2.0d0*oneWp**2+1.0d0)-2.0d0*oneWp**2*p*p
else
p = v1p(k)
discrim = (2.0d0*oneWp**2+1.0d0)-2.0d0*oneWp**2*p*p
endif
sqrtdiscrim = sqrt(discrim)
l_min_thick(2) = min((2.0d0*oneWp**2*p - sqrtdiscrim)/(2.0d0*oneWp**2+1.0d0),p)
l_max_thick(2) = max((2.0d0*oneWp**2*p + sqrtdiscrim)/(2.0d0*oneWp**2+1.0d0),p)
Expand All @@ -235,7 +259,11 @@ subroutine M1_explicitterms(dts,implicit_factor)
if (GR) then
discrim = alpp(k)**2*Xp(k)**2*3.0d0
else
discrim = 3.0d0
if (do_effectivepotential) then
discrim = alpp(k)**2*3.0d0
else
discrim = 3.0d0
endif
endif

sqrtdiscrim = sqrt(discrim)
Expand All @@ -254,7 +282,7 @@ subroutine M1_explicitterms(dts,implicit_factor)

!check for NaNs
if (l_min(1).ne.l_min(1)) then
write(*,*) "NaNs in speeds 1"
write(*,*) "NaNs in speeds 1",M1chi_space_minus(k+1),l_min_thin(1),l_min_thick(1),k,i,j
stop
else if(l_min(2).ne.l_min(2)) then
write(*,*) "NaNs in speeds 2"
Expand Down Expand Up @@ -301,10 +329,12 @@ subroutine M1_explicitterms(dts,implicit_factor)
endif

else
Jkplus1 = (1.0d0/(1.0d0-v1(k+1)**2))*(M1en_space(k+1)*(1.0d0+v1(k+1)**2* &
M1eddy_space(k+1))-2.0d0*M1flux_space(k+1)*v1(k+1))

Jk = (1.0d0/(1.0d0-v1(k)**2))*(M1en_space(k)*(1.0d0+v1(k)**2* &
M1eddy_space(k))-2.0d0*M1flux_space(k)*v1(k))
Jkplus1 = (1.0d0/(1.0d0-v1(k+1)**2))*(M1en_space(k+1)*(1.0d0+v1(k+1)**2* &
M1eddy_space(k+1))-2.0d0*M1flux_space(k+1)*v1(k+1))


if (a_asym.eq.1.0d0) then
diffusive_flux = 0.0d0
Expand Down Expand Up @@ -550,14 +580,28 @@ subroutine M1_explicitterms(dts,implicit_factor)
Xuprr(:)*X(k)**2*dWdr - Xuprr(:)*X(k)**4/alp(k)*dWvuprdt - &
heatterm_NL(:)/X(k)**2*dWvuprdr)
else
velocity_coeffs(:,1) = (oneW*(-2.0d0*Xupff(:)*onev*x1(k)) + &
Z(:)/alp(k)*dWdt + Yupr(:)*dWdr - &
Yupr(:)*dWvuprdt - Xuprr(:)*dWvuprdr)
if (do_effectivepotential) then
velocity_coeffs(:,1) = alp(k)*(oneW*(-2.0d0*Xupff(:)*onev*x1(k)) + &
Z(:)/alp(k)*dWdt + Yupr(:)*dWdr - &
Yupr(:)*dWvuprdt/alp(k) - Xuprr(:)*dWvuprdr)
velocity_coeffs(2,1) = velocity_coeffs(2,1) - alp(k)*dphidr(k)

!flux
velocity_coeffs(:,2) = alp(k)*(oneW*(-2.0d0*heattermff_NL(:)/x1(k)**2*onev*x1(k)) + &
Yupr(:)*dWdt/alp(k) + Xuprr(:)*dWdr - Xuprr(:)*dWvuprdt/alp(k) - &
heatterm_NL(:)*dWvuprdr)
velocity_coeffs(3,2) = velocity_coeffs(3,2) - alp(k)*dphidr(k)

!flux
velocity_coeffs(:,2) = (oneW*(-2.0d0*heattermff_NL(:)/x1(k)**2*onev*x1(k)) + &
Yupr(:)*dWdt + Xuprr(:)*dWdr - Xuprr(:)*dWvuprdt - &
heatterm_NL(:)*dWvuprdr)
else
velocity_coeffs(:,1) = (oneW*(-2.0d0*Xupff(:)*onev*x1(k)) + &
Z(:)/alp(k)*dWdt + Yupr(:)*dWdr - &
Yupr(:)*dWvuprdt - Xuprr(:)*dWvuprdr)

!flux
velocity_coeffs(:,2) = (oneW*(-2.0d0*heattermff_NL(:)/x1(k)**2*onev*x1(k)) + &
Yupr(:)*dWdt + Xuprr(:)*dWdr - Xuprr(:)*dWvuprdt - &
heatterm_NL(:)*dWvuprdr)
endif
endif

do i=1,number_species_to_evolve
Expand Down Expand Up @@ -728,14 +772,23 @@ subroutine M1_explicitterms(dts,implicit_factor)
X2 = X(k)**2
oneX = X(k)
else
invalp = 1.0d0
invalp2 = 1.0d0
alp2 = 1.0d0
onealp = 1.0d0
if (do_effectivepotential) then
invalp = 1.0d0/alp(k)
invalp2 = invalp**2
alp2 = alp(k)**2
onealp = alp(k)
else
invalp = 1.0d0
invalp2 = 1.0d0
alp2 = 1.0d0
onealp = 1.0d0
endif

invX = 1.0d0
invX2 = 1.0d0
X2 = 1.0d0
oneX = 1.0d0

endif

if (v_order.eq.-1) then
Expand Down
49 changes: 39 additions & 10 deletions src/M1/M1_implicitstep.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,14 +77,15 @@ subroutine M1_implicitstep(dts,implicit_factor)
integer :: info

logical :: nothappenyet1,nothappenyet2,stillneedconvergence
logical :: problem_fixing,trouble_brewing
logical :: problem_fixing,trouble_brewing,changedtwice
integer :: problem_zone
integer :: myloc(1)
real*8 :: maxRF

problem_fixing = .false.
problem_zone = 0
trouble_brewing = .true.
trouble_brewing = .false.
changedtwice = .false.

press_nu = 0.0d0
energy_nu = 0.0d0
Expand Down Expand Up @@ -265,7 +266,11 @@ subroutine M1_implicitstep(dts,implicit_factor)
sourceG(j,2) = implicit_factor*dts*onealp*W2* &
4.0d0*pi*x1(k)*rho(k)*h*(1.0d0+v2)
else
sourceG(j,2) = 0.0d0
if (do_effectivepotential) then
sourceG(j,2) = implicit_factor*dts*onealp*dphidr(k)
else
sourceG(j,2) = 0.0d0
endif
endif

!these are terms from the RHS of the flux equation
Expand Down Expand Up @@ -304,8 +309,13 @@ subroutine M1_implicitstep(dts,implicit_factor)
(press(k)+rho(k)*h*W2*v2)) - &
(eddyff(j)+eddytt(j))*invr)
else
sourceG(j+number_groups,1) = implicit_factor*dts*onealp*( -&
(eddyff(j)+eddytt(j))*invr)
if (do_effectivepotential) then
sourceG(j+number_groups,1) = implicit_factor*dts*onealp*( -&
(eddyff(j)+eddytt(j))*invr + dphidr(k))
else
sourceG(j+number_groups,1) = implicit_factor*dts*onealp*( -&
(eddyff(j)+eddytt(j))*invr)
endif
endif

enddo
Expand Down Expand Up @@ -1133,18 +1143,36 @@ subroutine M1_implicitstep(dts,implicit_factor)
nothappenyet2 = .false.
endif

if (changedtwice) then
stop "changed twice, must work unless explicit scattering is bad"
endif

!cases, low energy, energy coupling is large
!because difference between neighbouring bins is
!large, you can end up sending too much flux up in
!energy then is in the bin (hence, why we are
!here) Solution, make energy coupling term implicit.
if (j.le.3) then
sourceG(problem_zone,3) = q_M1_old(k,i,problem_zone,1) + &
B_M1(k,i,j,1) + D_M1(k,i,j,1)
sourceG(problem_zone,1) = sourceG(problem_zone,1) - &
C_M1(k,i,problem_zone,1)/q_M1_old(k,i,problem_zone,1)

!first time in here, try making energy coupling implcit
if (.not.trouble_brewing) then
sourceG(problem_zone,3) = q_M1_old(k,i,problem_zone,1) + &
B_M1(k,i,j,1) + D_M1(k,i,j,1)
sourceG(problem_zone,1) = sourceG(problem_zone,1) - &
C_M1(k,i,problem_zone,1)/q_M1_old(k,i,problem_zone,1)
endif

!second time in here, make flux implicit
if (trouble_brewing) then
sourceG(problem_zone,3) = q_M1_old(k,i,problem_zone,1) + &
D_M1(k,i,j,1)
sourceG(problem_zone,1) = sourceG(problem_zone,1) - &
B_M1(k,i,problem_zone,1)/q_M1_old(k,i,problem_zone,1)
changedtwice = .true.
endif

trouble_brewing = .true.
else if (j.ge.number_groups-3) then
else if (j.ge.number_groups-4) then
!sometimes the highest energy bin want to send
!out more energy flux than it has, fix that!
if (M1en_Exp_term(j).lt.0.0d0) then
Expand Down Expand Up @@ -1215,6 +1243,7 @@ subroutine M1_implicitstep(dts,implicit_factor)
enddo
problem_fixing = .false.
trouble_brewing = .false.
changedtwice = .false.
problem_zone = 0

!set and check new neutrino variables.
Expand Down
1 change: 1 addition & 0 deletions src/allocate_vars.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,7 @@ subroutine allocate_vars

allocate(ye(n1),ye_prev(n1))
allocate(dyedt_hydro(n1),dyedt_neutrino(n1))
allocate(depsdt(n1),dyedt(n1))
allocate(yep(n1))
allocate(yem(n1))
allocate(ynu(n1))
Expand Down
11 changes: 6 additions & 5 deletions src/eos.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -649,17 +649,18 @@ subroutine poly_eos(rho,enr,prs,dpdrho,dpde,soundsqr,keytemp)

! Local
!
prs = polyK * rho**polygamma
soundsqr = polygamma*polyK*rho**(polygamma-1.d0)

dpdrho = polyK * rho**(polygamma - 1) * polygamma
dpde = 0.0d0

if(keytemp.eq.1) then
! energy wanted
enr=prs/rho/(polygamma-1.d0)
soundsqr=polygamma*prs/rho
endif

prs = polyK * rho**polygamma
soundsqr = polygamma*polyK*rho**(polygamma-1.d0)

dpdrho = polyK * rho**(polygamma - 1) * polygamma
dpde = 0.0d0

end subroutine poly_eos

Expand Down
Loading

0 comments on commit dbffc85

Please sign in to comment.