Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

[FORTRAN] How did I screw up my code?

  1. Apr 15, 2015 #1

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    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:

    Code (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 helectron
    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.
     
  2. jcsd
  3. Apr 15, 2015 #2

    Mark44

    Staff: Mentor

    I count 567 lines of code.
    In the varelec_init() subroutine? Some detail on what "doesn't evolve correctly" means would be helpful. Presumably something untoward is happening after this code (if I'm in the right subroutine). In terms of the code and variables, what should be happening but isn't. You need to be more specific than the code "runs differently."
    Code (Fortran):
    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))
    .
    .
    .
    It's extremely difficult to follow this code as there is absolutely no indentation to make it easy to spot the bodies of each if.. elseif clause.
     
  4. Apr 15, 2015 #3

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    I can only see the end result, and that is the wave function. Basically what happens with varelecfrac equals 0 or 1 is that the terms in the wave function evolves slowly. When varelecfrac is set to 2 the wave function makes a sudden jump in evolution at around 901km (the simulation starts at 900km). The probability for a type 1 neutrino to evolve to a type 2 neutrino in the varelecfrac equals 0 or 1 case takes about 200-300km to grow from 0 to .001 whereas in the varelecfrac equals 2 case it takes 1km.

    The point is though that right now I have it set up so that no matter what I set varelecfrac (and spinco for that matter) equal, the code should be running (essentially) the same exact thing. But it is not.
     
  5. Apr 15, 2015 #4

    DrClaude

    User Avatar

    Staff: Mentor

    I'm a bit confused here. If I understand correctly, you mean this part of the code
    Code (Fortran):

      !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
     
    Shouldn't Yevar = 0.4d0 for all cases to be the same?

    By the way, the first part of elecfracsplint is very inefficient. Instead of going through the list of points until you find the where the desired position is with respect to the list, you should use a bisection method.
     
  6. Apr 15, 2015 #5

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    Well 1 or .4 it shouldn't make a huge difference, but yeah I set it equal to .4 as well and the code runs differently. I rewrote the code so it read like this:

    Code (Fortran):

      !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)  
      Yevar = .4d0
      endif
     
     
  7. Apr 15, 2015 #6

    DrClaude

    User Avatar

    Staff: Mentor

    But when varelecfrac = 1, Yevar is not necessarily 0.4.
     
  8. Apr 15, 2015 #7

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    Yeah it goes from .1 to .4...don't worry about the varelecfrac=1 for right now. The difference between .1 and .4 and 1 is significant, but not enough to make the wave function evolve 300+ times faster. varelecfrac=1 and varelecfrac=0 cases evolve similarly because the effects which are very sensitive to electron fraction I have not added to the code yet. The problem is that the code runs very very differently when varelecfrac=2, even if I manually set Yevar=.4.

    I first thought that I might have typo'd some of the arrays and perhaps overwrote my rhobgrid array with a yegrid array or something, but I looked through the code several times and did not find such an error.
     
  9. Apr 15, 2015 #8

    DrClaude

    User Avatar

    Staff: Mentor

    Then my suggestion is to comment out the call to elecfracsplint to see if anything nefarious is going on in that subroutine.
     
  10. Apr 15, 2015 #9

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    Ok, I commented out that line and the glitch is still present...:(

    My hunch is the error occurs in varelec_init because I've basically tried everything in the helectron and elecfracsplint to check to make sure it wasn't messing around with the actual electron fraction. I've even written code that basically just says "write(6,*) Yevar" inside helec and seen that the electron fraction is initially normal. The problem with this method is helec gets called tens of thousands of times at each radial step so if I put a write statement inside it and don't stop the code it will literally run me out of hard drive space saving this text file lol.

    EDIT: I checked Yevar at r=950 and it was still 1...so...
     
    Last edited: Apr 15, 2015
  11. Apr 15, 2015 #10

    Mark44

    Staff: Mentor

    Assuming that your hunch is correct about the error occurring in varelec_init(), it's not surprising to me that you get different behavior when varelecfrac is 2 as compared to when the values are 0 or 1. The only thing that is happening when varalecfrac is 0 or is 1 is that efname gets set to either 'Constant electron fraction' or 'Analytic profile', respectively. All the rest of the action in this subroutine occurs when varalecfrac is 2.

    Here is the block of code in varalec_init() starting with the check to see if varalecfrac equals 2.
    I have added indentation to your code to make it more readable.

    One thing that stands out to me is the backspace command 14 lines down. Per some Oracle docs for Fortran 77 that I found, "The BACKSPACE statement positions the specified file to just before the preceding record." Is there some reason for doing this? A comment would be useful giving an explanation for why this is needed.


    Code (Fortran):

      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
     
     
  12. Apr 15, 2015 #11

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    I'm not sure what the backspace is for, it was in the original code that I got... I just modified the portion to read also an electron fraction. I can comment that out and see what happens.

    EDIT: Ok, so I commented out the backspace and what happened was that the code read one less line of the electron fraction...I have 8 values in the grid but it read only the last 7.
     
    Last edited: Apr 15, 2015
  13. Apr 16, 2015 #12

    Mark44

    Staff: Mentor

    I'm not saying that the backspace thing was the root of your problem -- it's just that I didn't understand what its purpose was (and there are no comments to help a reader who is unfamiliar with the code).

    Did you understand the rest of my comment, the part where I said that virtually everything in varelec_init() is going on when varelecfrac is 2? The parts for varelecfrac values of 0 and 1 are essentially no-ops, doing nothing more than setting a couple of strings. It's no wonder that varelecfrac == 2 behaves differently than for the the other two values.
     
  14. Apr 16, 2015 #13

    DrClaude

    User Avatar

    Staff: Mentor

    Yes, but what varelec_init() does is set the values of ryegrid and yegrid for the interpolation of the electron fraction. If the interpolation is not used, as when elecfracsplint is commented out, there should be no difference. Anyway, a good way to check if the problem is in varelec_init is to change
    if(myid.eq.0) then
    to
    if(.false.) then
    such that that part of the code is not executed, in conjunction with commenting out the call to elecfracsplint.

    I notice also that the values of ryegrid and yegrid are not propagated to the other processors. Only myid = 0 will have the correct values.
     
  15. Apr 16, 2015 #14

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    You are right, I have not parallelized the code yet. I am currently running this code on a single core. I will comment out the call to varelec_init and see what happens.

    EDIT: Ok, so I commented out varelefrac_init and the call to elecfracsplint and set Yevar=1 and the bug went away. I think the problem is in varelecfrac_init then...not sure what though.. =/
     
    Last edited: Apr 16, 2015
  16. Apr 18, 2015 #15

    FactChecker

    User Avatar
    Science Advisor
    Gold Member

    A general debugging hint. Put write statements in at key points and examine the outputs till you can narrow the problem down. Concentrate on key interfaces to begin with and compare the working examples with the bad examples.
     
  17. Apr 18, 2015 #16

    mfb

    User Avatar
    2016 Award

    Staff: Mentor

    Another general hint: before you make signficant changes, make a backup (much better: use a version control system). If you really screw it up, you can look at all code changes, or go back to the last working version and implement new things in smaller steps to see what goes wrong.

    And can you fix Planck's name in variables please? Every future editor of that code will hate you for misspelled names as variables.

    The answer to the question "How did I screw up my code?" is in the title: FORTRAN
     
  18. Apr 19, 2015 #17

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    I never noticed that misspelling lol...but it's not mine, it was there when I got the code. I'm not sure which parts of the code uses that variable, so to fix this misspelling I'll have to look through all the other modules to figure out where I have to fix it. I'll get on that eventually.

    I think I have the problem narrowed down to the varelecfrac_init subroutine in the code, but as that subroutine should only initialize an electron fraction grid, I'm not sure what kind of write statements I can make to verify the variables. I have already written write statements for the electron fraction itself when the actual electron fraction is figured out in helec, and there are no problems with it, also the varelecfrac_init subroutine itself has write statements to show the grid that I inputted and there's no problems there...so I'm kind of at a loss. =/
     
  19. Apr 19, 2015 #18

    FactChecker

    User Avatar
    Science Advisor
    Gold Member

    @Matterwave , I owe you an apology. I see you already have a lot of write statements in your code. My comment was sloppy.
     
  20. Apr 20, 2015 #19

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    Oh my...so I actually edited out the part in varelecfrac_init that read in an electron fraction grid and just set the grid manually myself. The bug is exactly the same...

    I even cranked up the electron fraction to 100 and the results did not change. This is so baffling!

    It seems that if varelecfrac_init is called, no matter what I set the electron fraction to, the results are the same..
     
  21. Apr 20, 2015 #20

    Mark44

    Staff: Mentor

    There is no varelecfrac_init subroutine, at least in the code that you posted at the start of this thread. Is there such a sub in some other module that you didn't show us?
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook