Changes between Version 15 and Version 16 of u/bliu/pnfldiff

Timestamp:

08/10/18 16:34:48 (6 years ago)

Author:

Baowei Liu

Comment:

—

u/bliu/pnfldiff

@@ -20 +20 @@
 }}}
 
-== Basic Equations ==
+== Equations ==
 
 $ \frac{\partial n_{e}}{\partial t} + \nabla\cdot(n_{e} \bold{u}) = I-r $[[BR]]
@@ -60 +60 @@
 8. [https://astrobear.pas.rochester.edu/trac/wiki/u/adebrech/Matlab/RadiationPressure Radiation pressure 2]
 
-8. [https://astrobear.pas.rochester.edu/trac/wiki/u/adebrech/Matlab/IonizationFront Ionization front]
+9. [https://astrobear.pas.rochester.edu/trac/wiki/u/adebrech/Matlab/IonizationFront Ionization front]
+
+== Current code ==
+
+{{{
+! Goal is to take input flux and integrate equations across grid to find ionization rate and heating rate
+  ! Flux has units of number per area per time in computational units 
+  ! Data has computational units of the density of neutral hydrogen, ionized hydrogen, ionization rate, and heating rate
+  SUBROUTINE DoLineTransfer(data, rates)
+    INTEGER :: i,j,k
+    REAL(KIND=qPREC) :: emtau, f, dx, tau, dtau
+    REAL(KIND=qpREC), DIMENSION(:,:,:,:) :: data, rates
+    !Scale parameters to computational scales
+
+    dx=levels(MaxLevel)%dx
+!    write(*,*) LineTransferFlux
+
+    !$OMP PARALLEL DO SCHEDULE(STATIC) COLLAPSE(2) SHARED(sigmaH_scale, egamma_scale, dx, rates, data, LineTransferFlux) PRIVATE(i,j,k,tau,f,dtau,emtau) DEFAULT(NONE)
+    DO j=1,size(data,2)
+       DO k=1,size(data,3)
+          tau=0d0
+          f=LineTransferFlux
+         !write(*,'(A4,20A25)'), 'i', 'f', 'nH','nHII', 'exp(-tau)', 'dtau', 'sigmaH_scale', 'dx'
+          DO i=1,size(data,1)
+
+             !First calculate exp(-tau) where tau =sigmaH*nH+sigmad*nHII
+!             emtau=exp(-(sigmaH_scale*data(i,j,k,iH) + sigmad_scale*data(i,j,k,iHII))*dx) )
+             dtau=(sigmaH_scale*data(i,j,k,lt_iH)*dx)
+             emtau=exp(-dtau)
+             !if (emtau /= 1.) write(*,'(I4,20E25.15)') i, f, data(i,j,k,iH:iHII), emtau,(sigmaH_scale*data(i,j,k,iH))*dx, sigmaH_scale, dx
+             !Store the photon absorption rate used to calculate the heating and ionization rates
+             rates(i,j,k,iioniz)=f*(1d0-emtau)/dx               !has computational units of number per volume per second
+!             rates(i,j,k,iheat)=0
+             rates(i,j,k,iheat)=f*(1d0-emtau)/dx*egamma_scale  !has computational units of energy per volume per second
+
+             !Update the flux for the next cell along this ray
+             f=f*emtau
+             tau=tau+dtau
+             data(i,j,k,lt_iD)=tau
+          END DO
+       END DO
+    END DO
+    !$OMP END PARALLEL DO
+  END SUBROUTINE DoLineTransfer
+}}}
+
+