- #1
- 3,971
- 328
Hi guys, there's a bug in my code which I've been trying to get rid of for the last 2 weeks. I'm completely out of ideas now on what it is so I have to turn to you guys for help. The code is here:
Sorry for the long piece of code.
Basically what I've added to this piece of code is the ability to have a variable electron fraction. The integer varelecfrac is an input which can be 0, 1, or 2 depending on what kind of a variable electron fraction I want to have. My code runs fine if I set varelecfrac equal to 0 or 1, but when I set it equal to 2 it doesn't evolve correctly.
I am completely baffled because I set the electron fraction grid to be 1 for all radii but the code still runs differently if I set varelecfrac equal to 2 than equal to 1 or 0 (in both scenarios the electron fraction has been set to 1 for testing purposes).
In fact, I even wrote a line Yevar=1.d0 right after I call elecfracsplint, thereby overriding all ability for the code to mess around with Yevar and the code STILL runs differently.
I am totally baffled at this point. I don't even know what else to try debugging.
iprofile is set to 3, spinco is set to 0 or 1 and the problem persists regardless of what I set spinco to. If I set spinco to 0, the code should actually not worry about this variable electron fraction at all since in that case it should be reading an electron density grid directly instead of a baryon density grid. But this problem persists! What is happening...T_T
I apologize for the implicit variables in this code. It was written this way when I got it.
Fortran:
module helectron
use params
use parallel
!* electron density parameters
integer :: iprofile
!* Ref: Full & Qian, astro-ph/0505240, Eqs. (71-74)
! solar mass in kg
real(kind=reel8),parameter :: Msun = 1.98844d+30
! proton mass in kg
real(kind=reel8),parameter :: Mproton = 1.67262171d-27
! plank mass in kg
real(kind=reel8),parameter :: Mplank = 2.17645d-8
! mass of the neutron star in kg
real(kind=reel8),parameter :: Mns = 1.4*Msun
! entropy per baryon in k_B
real(kind=reel8),parameter :: Sb = 1.4d+2
! radius of neutrinosphere in cm
real(kind=reel8),parameter :: Rnu = 1.1d+6
! temperature in cm**-1 at neutrinosphere
real(kind=reel8),parameter :: T0cm = (Mproton*Mns)/(Mplank*Mplank) &
& /Sb/Rnu
! electron fraction, Yevar used for variable electron fraction
real(kind=reel8),parameter :: Ye = 1.d0
real(kind=reel8) :: Yevar
! electron density at neutrino sphere in cm**-3
real(kind=reel8) :: ne0 = Ye*(2.d0*pi*pi/45.d0)*T0cm**3/Sb
! baryon density at neutrino sphere in cm**-3
real(kind=reel8) :: nb0 = (2.d0*pi*pi/45.d0)*T0cm**3/Sb
!
! grid for electron density, baryon density and electron fraction
integer :: iegrid
real(kind=reel8),dimension(1000) :: rgrid,ryegrid,rhoegrid,rhoe2grid
real(kind=reel8),dimension(1000) :: rhobgrid,rhob2grid,yegrid,ye2grid
!
! mixing angles and masses
!
real(kind=reel8) :: dm21sq,dm31sq,dm32sq
real(kind=reel8) :: theta12, theta13, theta23,deltacp
real(kind=reel8),dimension(nflavor,nflavor) :: hvac,hmix,hmass
contains
subroutine hvac_init
if(nflavor.eq.2) then
dm20= dm21sq
cos2v = cos( 2.d0*theta13)
sin2v = sin( 2.d0*theta13)
if(myid.eq.0) then
write(6,'(a20,1pe15.7) ') '2-flavor vacuum matrix'
write(6,'(a20,1pe15.7) ') 'dm20 = ',dm20
write(6,'(a20,1pe15.7) ') 'cos2v = ',cos2v
write(6,'(a20,1pe15.7) ') 'sin2v = ',sin2v
endif
elseif(nflavor.eq.3) then
c13=cos(theta13)
s13=sin(theta13)
c12=cos(theta12)
s12=sin(theta12)
c23=cos(theta23)
s23=sin(theta23)
hmix(1,1)=c12*c13
hmix(1,2)=s12*c13
hmix(1,3)=s13
hmix(2,1)=-s12*c23-c12*s23*s13
hmix(2,2)=c12*c23-s12*s23*s13
hmix(2,3)=s23*c13
hmix(3,1)=s12*s23-c12*c23*s13
hmix(3,2)=-c12*s23-s12*c23*s13
hmix(3,3)=c23*c13
hmass(:,:)=0.d0
hmass(2,2)=dm21sq
hmass(3,3)=dm21sq+dm32sq
hvac(:,:)=0.d0
do i=1,nflavor
do j=1,nflavor
do k=1,nflavor
do l=1,nflavor
hvac(i,l)=hvac(i,l)+hmix(i,j)*hmass(j,k)*hmix(l,k)
enddo
enddo
enddo
enddo
if(myid.eq.0) then
write(6,'(a20,1pe15.7) ') '3-flavor vacuum matrix'
write(6,'(a20,1pe15.7) ') 'dm21sq = ',dm21sq
write(6,'(a20,1pe15.7) ') 'dm32sq = ',dm32sq
write(6,'(a20,1pe15.7) ') 'theta12 = ',theta12
write(6,'(a20,1pe15.7) ') 'theta13 = ',theta13
write(6,'(a20,1pe15.7) ') 'theta23 = ',theta23
write(6,'(a20,1pe15.7) ') 'cos(theta12) = ',c12
write(6,'(a20,1pe15.7) ') 'sin(theta12) = ',s12
write(6,'(a20,1pe15.7) ') 'cos(theta13) = ',c13
write(6,'(a20,1pe15.7) ') 'sin(theta13) = ',s13
write(6,'(a20,1pe15.7) ') 'cos(theta23) = ',c23
write(6,'(a20,1pe15.7) ') 'sin(theta23) = ',s23
write(6,'(a20,1pe15.7) ') 'deltacp = ',0.
write(6,'(a40,1pe15.7) ') '3-flavor mixing matrix '
write(6,'(1p,3e15.7) ') hmix(1,1),hmix(1,2),hmix(1,3)
write(6,'(1p,3e15.7) ') hmix(2,1),hmix(2,2),hmix(2,3)
write(6,'(1p,3e15.7) ') hmix(3,1),hmix(3,2),hmix(3,3)
endif
! Do some unitarity checks on mixing matrix
sumsq1 = sum(hmix(1,1:3)**2)-1.d0
if(abs(sumsq1).gt.1.d-5) then
write(6,*) ' Error sumsq1 = ', sumsq1
stop
endif
sumsq2 = sum(hmix(2,1:3)**2)-1.d0
if(abs(sumsq2).gt.1.d-5) then
write(6,*) ' Error sumsq2 = ', sumsq2
stop
endif
sumsq3 = sum(hmix(3,1:3)**2)-1.d0
if(abs(sumsq3).gt.1.d-5) then
write(6,*) ' Error sumsq3 = ', sumsq3
stop
endif
traceo3=(hvac(1,1)+hvac(2,2)+hvac(3,3))/3.d0
hvac(1,1)=hvac(1,1)-traceo3
hvac(2,2)=hvac(2,2)-traceo3
hvac(3,3)=hvac(3,3)-traceo3
if(myid.eq.0) then
write(6,'(a40,1pe15.7) ') '3-flavor vacuum Hamiltonian (traceless)'
write(6,'(1p3e15.7) ') hvac(1,1),hvac(1,2),hvac(1,3)
write(6,'(1p3e15.7) ') hvac(2,1),hvac(2,2),hvac(2,3)
write(6,'(1p3e15.7) ') hvac(3,1),hvac(3,2),hvac(3,3)
endif
endif
return
end subroutine hvac_init
subroutine helec_init(iprofile,spinco,filename,helecname)
character*80 filename
character*80 helecname
character*80 rest
character*1 first
integer :: spinco
#ifdef mpi
include 'mpif.h'
#endif
helecname = ''
! iprofile=1 or iprofile=2 are simple analytic profiles
if(iprofile.eq.1) then
helecname = 'Fuller and Qian 1/r**3'
elseif(iprofile.eq.2) then
helecname = 'Fuller and Qian 1/r**3 + exp '
elseif(iprofile.eq.3) then
if(spinco.eq.0) then
if(myid.eq.0) then
write(6,*) ' Reading electron density from ',filename
write(6,*) ' Beginning of file '
icomments=1
open(unit=20,file=trim(filename))
do while (icomments.eq.1)
read(20,'(a1,a80)')first,rest
if(helecname.eq.'') helecname=rest
if (first.eq.'#') then
write(6,'(a1,a80)') first,rest
else
icomments=0
backspace 20
endif
enddo
! read data file
iend=0
iegrid=0
do while (iend.eq.0)
iegrid=iegrid+1
read(20,*,iostat=ierr) rgrid(iegrid), rhoegrid(iegrid)
if(ierr.ne.0) then
iend=1
iegrid=iegrid-1
close(20)
endif
enddo
write(6,*) ' Electron density grid contains ',iegrid,'values'
write(6,*) ' Beginning of rho_e grid '
do i=1,min(5,iegrid)
write(6,'(f15.7,e18.9)') rgrid(i),rhoegrid(i)
enddo
write(6,*) ' End of rho_e grid '
do i=max(1,iegrid-4),iegrid
write(6,'(f15.7,e18.9)') rgrid(i),rhoegrid(i)
enddo
! convert rhoegrid from ln(n_e*km^3) to ln (ne_e*cm^3)
! conversion factor is 10**(-15)
write(6,*) ' initial density = ',exp(rhoegrid(1)),' / km**3 '
rhoegrid(1:iegrid)=rhoegrid(1:iegrid)-log(1.d15)
write(6,*) ' converting to ln(rho_e*cm**3) = ',exp(rhoegrid(1)),' /cm**3'
! end of if for myid = 0
endif
#ifdef mpi
call mpi_bcast(iegrid,1,mpi_integer,0,mpi_comm_world,ierr)
call mpi_bcast(rgrid,iegrid,mpi_double_precision,0,mpi_comm_world,ierr)
call mpi_bcast(rhoegrid,iegrid,mpi_double_precision,0,mpi_comm_world,ierr)
#endif
if(myid.eq.0) then
write(6,*) ' setting up spline'
endif
call spline(rgrid,rhoegrid,iegrid,1.d31,1.d31,rhoe2grid)
if(myid.eq.0) then
write(6,*) ' done setting up spline'
endif
elseif(spinco.eq.1.or.spinco.eq.2) then
if(myid.eq.0) then
write(6,*) ' Reading baryon density from ',filename
write(6,*) ' Beginning of file '
icomments=1
open(unit=20,file=trim(filename))
do while (icomments.eq.1)
read(20,'(a1,a80)')first,rest
if(helecname.eq.'') helecname=rest
if (first.eq.'#') then
write(6,'(a1,a80)') first,rest
else
icomments=0
backspace 20
endif
enddo
! read data file
iend=0
iegrid=0
do while (iend.eq.0)
iegrid=iegrid+1
read(20,*,iostat=ierr) rgrid(iegrid), rhobgrid(iegrid)
if(ierr.ne.0) then
iend=1
iegrid=iegrid-1
close(20)
endif
enddo
write(6,*) ' Baryon density grid contains ',iegrid,'values'
write(6,*) ' Beginning of rho_b grid '
do i=1,min(5,iegrid)
write(6,'(f15.7,e18.9)') rgrid(i),rhobgrid(i)
enddo
write(6,*) ' End of rho_b grid '
do i=max(1,iegrid-4),iegrid
write(6,'(f15.7,e18.9)') rgrid(i),rhobgrid(i)
enddo
! convert rhobgrid from ln(n_b*km^3) to ln (nb_e*cm^3)
! conversion factor is 10**(-15)
write(6,*) ' initial density = ',exp(rhobgrid(1)),' / km**3 '
rhobgrid(1:iegrid)=rhobgrid(1:iegrid)-log(1.d15)
write(6,*) ' converting to ln(rho_b*cm**3) = ',exp(rhobgrid(1)),' /cm**3'
! end of if for myid = 0
endif
#ifdef mpi
call mpi_bcast(iegrid,1,mpi_integer,0,mpi_comm_world,ierr)
call mpi_bcast(rgrid,iegrid,mpi_double_precision,0,mpi_comm_world,ierr)
call mpi_bcast(rhoegrid,iegrid,mpi_double_precision,0,mpi_comm_world,ierr)
#endif
if(myid.eq.0) then
write(6,*) ' setting up spline'
endif
call spline(rgrid,rhobgrid,iegrid,1.d31,1.d31,rhob2grid)
if(myid.eq.0) then
write(6,*) ' done setting up spline'
endif
! end of if for spinco = 1 or 2
endif
! end of if for iprofile = 1,2, or 3
endif
return
end subroutine helec_init subroutine varelec_init(varelecfracf,elecprofilef,filename,efname)
integer,intent(in) :: varelecfracf,elecprofilef
character*80 filename
character*80 efname
character*80 rest
character*1 first efname = ''
! varelecfracf = 0 means constant elecfrac, 1 means analytic profile, 2 means read from file
! elecprofilef = 1 means linear electron fraction
if(varelecfracf.eq.0) then
efname = 'Constant electron fraction'
elseif(varelecfracf.eq.1) then
efname = 'Analytic profile '
elseif(varelecfracf.eq.2) then
if(myid.eq.0) then
write(6,*) ' Reading electron fraction from ',filename
write(6,*) ' Beginning of file '
icomments=1
open(unit=20,file=trim(filename))
do while (icomments.eq.1)
read(20,'(a1,a80)')first,rest
if(efname.eq.'') efname=rest
if (first.eq.'#') then
write(6,'(a1,a80)') first,rest
else
icomments=0
backspace 20
endif
enddo
! read data file
iend=0
iegrid=0
do while (iend.eq.0)
iegrid=iegrid+1
read(20,*,iostat=ierr) ryegrid(iegrid), yegrid(iegrid)
if(ierr.ne.0) then
iend=1
iegrid=iegrid-1
close(20)
endif
enddo
write(6,*) ' Electron fraction grid contains ',iegrid,'values'
write(6,*) ' Beginning of Y_e grid '
do i=1,min(5,iegrid)
write(6,'(f15.7,e18.9)') ryegrid(i),yegrid(i)
enddo
write(6,*) ' End of Y_e grid '
do i=max(1,iegrid-4),iegrid
write(6,'(f15.7,e18.9)') ryegrid(i),yegrid(i)
enddo
! end of if for myid = 0
endif! end of if for varelecfracf=0,1,2
endif
return
end subroutine varelec_init
subroutine elecfracsplint(ryegrid,yegrid,r,Yevar)
real(kind=reel8),intent(in) :: r
real(kind=reel8),intent(out) :: Yevar
real(kind=reel8),dimension(1000),intent(in) :: ryegrid,yegrid
integer :: igrid
igrid=0
do
if(r.ge.ryegrid(igrid)) then
igrid=igrid+1
if(igrid.ge.999) then
write(6,*) ' Error: Exceeded electron fraction grid '
stop
endif
else
exit
endif
enddo
Yevar=(yegrid(igrid)-yegrid(igrid-1))*r/(ryegrid(igrid)-ryegrid(igrid-1))+yegrid(igrid)- &
& (ryegrid(igrid)*yegrid(igrid)-ryegrid(igrid)*yegrid(igrid-1))/ &
& (ryegrid(igrid)-ryegrid(igrid-1))
return
end subroutine function helec(r,e,iphi,iprofile,spinco,varelecfrac)
implicit real(kind=reel8) (a-h,o-z)
real(kind=reel8),dimension(nflavor,nflavor,2) :: helec
integer :: spinco,varelecfrac,iprofile
real(kind=reel8) :: he,hence
real(kind=reel8) :: Yevar
! henc is neutral current term that must be added back in if spin coherence is used
real(kind=reel8),intent(IN) :: r,e
real(kind=reel8) :: gs !statistical weight
real(kind=reel8) :: r3 ! (r/Rnu)**3
! integer init/0/
! save init
helec(:,:,:)=(0.d0,0.d0)
Yevar=.4
!initialize variable electron fraction
if(varelecfrac.eq.0) then
Yevar = .4d0
elseif(varelecfrac.eq.1) then
Yevar = elecfrac(r,1)
elseif(varelecfrac.eq.2) then
call elecfracsplint(ryegrid,yegrid,r,Yevar)
endif!* Ref: Full & Qian, astro-ph/0505240, Eqs. (71-74)
if(iprofile.eq.1.or.iprofile.eq.2) then
gs = 5.5d0 !statistical weight
r3 = (r*1.d5/Rnu)**3
if(varelecfrac.eq.0) then
he = 2.d0*gs*gfermiort2*ne0/r3
elseif(varelecfrac.eq.1.or.varelecfrac.eq.2) then
he = 2.d0*gs*gfermiort2*nb0*Yevar/r3
endif
if(spinco.eq.1.or.spinco.eq.2) then
henc = gs*gfermiort2*nb0*(1.d0-Yevar)/r3
endif
endif
! for iprofile=2 add an exponential decaying term
! with large rho_e near surface
if(iprofile.eq.2) then
!--------------------------
! electron density in g/cm**3
heexp=2.7d12*exp(- (r-11.d0)/0.18d0)
! multiply by
! electron fraction = 0.4
! Avagadro's number= 6.023d23=ava
! to get density in 1/cm**3
if(varelecfrac.eq.0) then
he=2.d0*heexp*0.4d0*ava*gfermiort2 + he
elseif(varelecfrac.eq.1.or.varelecfrac.eq.2) then
he=2.d0*heexp*Yevar*ava*gfermiort2 + he
endif
if(spinco.eq.1.or.spinco.eq.2) then
henc = heexp*ava*gfermiort2*(1.d0-Yevar) + henc
endif
endif!--------------------------
if(iprofile.eq.3) then
rlog=log(r)
if(spinco.eq.0) then
call splint(rgrid,rhoegrid,rhoe2grid,iegrid,rlog,rhoelog)
he=2.d0*gfermiort2*exp(rhoelog)
elseif(spinco.eq.1.or.spinco.eq.2) then
call splint(rgrid,rhobgrid,rhob2grid,iegrid,rlog,rhoblog)
if(varelecfrac.eq.0) then
he=2.d0*gfermiort2*exp(rhoblog)*Ye
henc=gfermiort2*exp(rhoblog)*(1.d0-Ye)
elseif(varelecfrac.eq.1.or.varelecfrac.eq.2) then
he=2.d0*gfermiort2*exp(rhoblog)*Yevar
henc=gfermiort2*exp(rhoblog)*(1.d0-Yevar)
endif
endif
endif
if(nflavor.eq.2) then
hm=dm20/(4.d0*e)*cos2v
hm2=dm20/(4.d0*e)*sin2v
if(spinco.eq.0) then
helec(1,1,1)=.5d0*he-hm
helec(1,2,1)=hm2
helec(2,1,1)=hm2
helec(2,2,1)=-.5d0*he+hm
helec(1,1,2)=-.5d0*he-hm
helec(1,2,2)=hm2
helec(2,1,2)=hm2
helec(2,2,2)=.5d0*he+hm
! the factors of .5 come from removing the trace
elseif(spinco.eq.1) then
helec(1,1,1)=he-hm-henc
helec(1,2,1)=hm2
helec(2,1,1)=hm2
helec(2,2,1)=hm-henc
helec(1,1,2)=-he-hm+henc
helec(1,2,2)=hm2
helec(2,1,2)=hm2
helec(2,2,2)=hm+henc
! can not "remove trace" as the spin coherence Hamiltonian is already traceless
endif elseif(nflavor.eq.3) then
! 1 is electron flavor neutrino
! 2 is mu flavor neutrino
! 3 is tau flavor neutrino
helec(:,:,1)=hvac(:,:)/(2.d0*e)
helec(:,:,2)=hvac(:,:)/(2.d0*e)
if(spinco.eq.0) then
! add electron potential, keep it traceless
helec(1,1,1)=helec(1,1,1)+2.*he/3.d0
helec(2,2,1)=helec(2,2,1)- he/3.d0
helec(3,3,1)=helec(3,3,1)- he/3.d0
helec(1,1,2)=helec(1,1,2)-2.*he/3.d0
helec(2,2,2)=helec(2,2,2)+ he/3.d0
helec(3,3,2)=helec(3,3,2)+ he/3.d0
!
! if(init.eq.0) then
! write(6,*) ' Hamiltonian, R,E,he = ',r,e,he
! write(6,*) ' Neutrinos '
! write(6,'(1p,3e17.9)') helec(1,1:3,1)
! write(6,'(1p,3e17.9)') helec(2,1:3,1)
! write(6,'(1p,3e17.9)') helec(3,1:3,1)
! write(6,*) ' Anti-Neutrinos '
! write(6,'(1p,3e17.9)') helec(1,1:3,2)
! write(6,'(1p,3e17.9)') helec(2,1:3,2)
! write(6,'(1p,3e17.9)') helec(3,1:3,2)
! init=1
! endif
elseif(spinco.eq.1) then
! again, this hamiltonian will already be traceless
helec(1,1,1)=helec(1,1,1)+he-henc
helec(2,2,1)=helec(2,2,1)-henc
helec(3,3,1)=helec(3,3,1)-henc
helec(1,1,2)=helec(1,1,2)-he+henc
helec(2,2,2)=helec(2,2,2)+henc
helec(3,3,2)=helec(3,3,2)+henc
endif endif
return
end function helec
function elecfrac(r,elecprofile)
! We can choose from a variety of functional forms for the electron fraction
implicit real(kind=reel8) (a-h,o-z)
real(kind=reel8) :: elecfrac
real(kind=reel8),intent(in) :: r
integer,intent(in) :: elecprofile
if(elecprofile.eq.1) then
if(r.lt. 9.d+2) then
elecfrac=.1d0
elseif(r.ge. 9.d+2 .and. r.le. 5.d+3) then
elecfrac=(.3d0*r/4.1d+3)+.03415d0
elseif(r.gt. 5.d+3) then
elecfrac=.4d0
endif
elseif(elecprofile.eq.2) then
elecfrac=.4d0
else
elecfrac=.4d0
endif
return
end function elecfrac
end module helectronSorry for the long piece of code.
Basically what I've added to this piece of code is the ability to have a variable electron fraction. The integer varelecfrac is an input which can be 0, 1, or 2 depending on what kind of a variable electron fraction I want to have. My code runs fine if I set varelecfrac equal to 0 or 1, but when I set it equal to 2 it doesn't evolve correctly.
I am completely baffled because I set the electron fraction grid to be 1 for all radii but the code still runs differently if I set varelecfrac equal to 2 than equal to 1 or 0 (in both scenarios the electron fraction has been set to 1 for testing purposes).
In fact, I even wrote a line Yevar=1.d0 right after I call elecfracsplint, thereby overriding all ability for the code to mess around with Yevar and the code STILL runs differently.
I am totally baffled at this point. I don't even know what else to try debugging.
iprofile is set to 3, spinco is set to 0 or 1 and the problem persists regardless of what I set spinco to. If I set spinco to 0, the code should actually not worry about this variable electron fraction at all since in that case it should be reading an electron density grid directly instead of a baryon density grid. But this problem persists! What is happening...T_T
I apologize for the implicit variables in this code. It was written this way when I got it.