Merge pull request #7 from evanoconnor/mom_nu

including neutrino momentum into effective potential
evanoconnor · Jul 13, 2017 · 1e20db1 · 1e20db1
2 parents 5376470 + f89c8ad
commit 1e20db1
Show file tree

Hide file tree

Showing 6 changed files with 24 additions and 5 deletions.
diff --git a/src/GR1D_module.F90 b/src/GR1D_module.F90
@@ -221,6 +221,7 @@ module GR1D_module
   ! specific internal energy, cell center and at interfaces
   real*8,allocatable,save :: eps(:),epsp(:),epsm(:)
   real*8,allocatable,save :: energy_nu(:)
+  real*8,allocatable,save :: mom_nu(:)
   ! mass interior cell center, mass of cell, volume of cell
   real*8,allocatable,save :: mass(:), mass1(:), volume(:)
   ! pressure 

diff --git a/src/M1/M1_implicitstep.F90 b/src/M1/M1_implicitstep.F90
@@ -89,6 +89,7 @@ subroutine M1_implicitstep(dts,implicit_factor)
 
   press_nu = 0.0d0
   energy_nu = 0.0d0
+  mom_nu = 0.0d0
   ynu = 0.0d0
 
   nothappenyet1 = .true.
@@ -1338,6 +1339,7 @@ subroutine M1_implicitstep(dts,implicit_factor)
 
            press_nu(k) = press_nu(k) + oneeddy*oneM1en*4.0d0*pi*invX2**2
            energy_nu(k) = energy_nu(k) + oneM1en*4.0d0*pi
+           mom_nu(k) = mom_nu(k) + oneM1flux*4.0d0*pi
            !$OMP END CRITICAL        
 
         enddo

diff --git a/src/allocate_vars.F90 b/src/allocate_vars.F90
@@ -35,6 +35,7 @@ subroutine allocate_vars
   allocate(epsm(n1))
 
   allocate(energy_nu(n1))
+  allocate(mom_nu(n1))
 
   allocate(eps_kin(n1))
   allocate(binding_energy(n1))

diff --git a/src/initialize_vars.F90 b/src/initialize_vars.F90
@@ -133,6 +133,7 @@ subroutine initialize_arrays
   epsm(:) = 0.0d0
 
   energy_nu(:) = 0.0d0
+  mom_nu(:) = 0.0d0
 
   eps_kin(:) = 0.0d0
   binding_energy(:) = 0.0d0

diff --git a/src/mass_interior.F90 b/src/mass_interior.F90
@@ -43,7 +43,7 @@ subroutine mass_interior
   else
      if (do_effectivepotential) then
         mass(ghosts1+1) = 4.0d0*pi/3.0d0*x1(ghosts1+1)**3*(rho(ghosts1+1)* &
-             (1.0d0 + eps(ghosts1+1)) + energy_nu(ghosts1+1)*include_nus)
+             (1.0d0 + eps(ghosts1+1)) + (energy_nu(ghosts1+1)+mom_nu(ghosts1+1)*v1(ghosts1+1))*include_nus)
         mass(ghosts1+1) = mass(ghosts1+1)*sqrt(1.0d0-2.0d0*mass(ghosts1+1)/x1(ghosts1+1))
 
         dphidr(ghosts1+1) = (mass(ghosts1+1) + 4.0d0*pi*x1(ghosts1+1)**3* &
@@ -53,21 +53,21 @@ subroutine mass_interior
 
         do i=ghosts1+2,n1-1
            mass(i) = mass(i-1) + &
-                4.0d0/3.0d0*pi*(rho(i-1)*(1.0d0+eps(i-1))+energy_nu(i-1)*include_nus) &
+                4.0d0/3.0d0*pi*(rho(i-1)*(1.0d0+eps(i-1))+(energy_nu(i-1)+mom_nu(i-1)*v1(i-1))*include_nus) &
                 * ( x1i(i)**3 - x1(i-1)**3 )*sqrt(1.0d0-2.0d0*mass(i-1)/x1(i-1))
            mass(i) = mass(i) + &
-                4.0d0/3.0d0*pi*(rho(i)*(1.0d0+eps(i))+energy_nu(i)*include_nus) * &
+                4.0d0/3.0d0*pi*(rho(i)*(1.0d0+eps(i))+(energy_nu(i)+mom_nu(i)*v1(i))*include_nus) * &
                 (x1(i)**3 - x1i(i)**3)*sqrt(1.0d0-2.0d0*mass(i)/x1i(i))
            dphidr(i) = (mass(i) + 4.0d0*pi*x1(i)**3*(press(i)+press_nu(i)))/ &
                 (x1(i)**2*(1.0d0+v1(i)**2-2.0d0*mass(i)/x1(i)))* &
                 (rho(i)+eps(i)*rho(i)+press(i))/rho(i)
         enddo
 
         mass(n1) = mass(n1-1) + &
-             4.0d0/3.0d0*pi*(rho(n1-1)*(1.0d0+eps(n1-1))+energy_nu(n1-1)*include_nus)* &
+             4.0d0/3.0d0*pi*(rho(n1-1)*(1.0d0+eps(n1-1))+(energy_nu(n1-1)+mom_nu(n1-1)*v1(n1-1))*include_nus)* &
              (x1i(n1)**3 - x1(n1-1)**3)*sqrt(1.0d0-2.0d0*mass(n1-1)/x1(n1-1))
         mass(n1) = mass(n1) + &
-             4.0d0/3.0d0*pi*(rho(n1)*(1.0d0+eps(n1))+energy_nu(n1)*include_nus)* &
+             4.0d0/3.0d0*pi*(rho(n1)*(1.0d0+eps(n1))+(energy_nu(n1)+mom_nu(n1)*v1(n1))*include_nus)* &
              (x1(n1)**3 - x1i(n1)**3)*sqrt(1.0d0-2.0d0*mass(n1)/x1i(n1))         
         dphidr(n1) = (mass(n1) + 4.0d0*pi*x1(n1)**3*(press(n1)+press_nu(n1)))/ &
              (x1(n1)**2*(1.0d0+v1(n1)**2-2.0d0*mass(n1)/x1(n1)))* &

diff --git a/src/restart_H5.F90 b/src/restart_H5.F90
@@ -548,6 +548,15 @@ subroutine restart_output_h5
      call h5dclose_f(dset_id, error)
      call h5sclose_f(dspace_id, error)
      cerror = cerror + error
+
+     call h5screate_simple_f(rank, dims1, dspace_id, error)
+     call h5dcreate_f(file_id, "mom_nu", H5T_NATIVE_DOUBLE,&
+          & dspace_id, dset_id, error)
+     call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, mom_nu, dims1,&
+          & error)
+     call h5dclose_f(dset_id, error)
+     call h5sclose_f(dspace_id, error)
+     cerror = cerror + error
 
      call h5screate_simple_f(rank, dims1, dspace_id, error)
      call h5dcreate_f(file_id, "dnupdr", H5T_NATIVE_DOUBLE, dspace_id&
@@ -1076,6 +1085,11 @@ subroutine restart_init_h5
      call h5dclose_f(dset_id,error) 
      cerror = cerror + error
 
+     call h5dopen_f(file_id, "mom_nu", dset_id, error)
+     call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, mom_nu, dims1, error)
+     call h5dclose_f(dset_id,error) 
+     cerror = cerror + error
+
      call h5dopen_f(file_id, "dnupdr", dset_id, error)
      call h5dread_f(dset_id, H5T_NATIVE_DOUBLE, dnupdr, dims1, error)
      call h5dclose_f(dset_id,error)