Project

Like to use the Flux Lim Diffusion to track ionization fronts
I know this won’t work well but its better than nothing.

Jonathan what would it take to do a test of this?

Ionizing source in media with initially low optical depth but tau builds up till it hits 1
We are basically trying to do the strongmen(Stromgren?) sphere problem.

tau → 1, exp(-tau) or ionization ratio??

  I think the current FLD Radiation transport (or the thermal diffusion) scheme could be adapted.  You would have to sub-cycle within the ionization/ radiation solve as the ionization fraction changes - though we do that already with the other parabolic solvers.  Usually there's some additional restriction on the hydro time step - that we would need to work out - or find in the literature...

But a first step would be to gather the necessary equations, linearize them, and then discretize them...  similar to what is done for the thermal diffusion solver.  Using the same radiation flux limiter would be straight-forward.  

https://astrobear.pas.rochester.edu/trac/wiki/ThermalConduction

Equations

\frac{\partial n_{e}}{\partial t} + \nabla\cdot(n_{e} \bold{u}) = I-r

I is the rate of ionizations per unit volume and r is the rate of recombinations per unit volume

I = -\nabla\cdot J

J is the flux of ionizing photons. Inside the front where the gas is nearly completely ionized.

The recombination rate

r=\alpha n_{e}^{2}

at T~10000K, \alpha=2.6\times10^{-13} cm^{3}s^{-1}

Optical depth of the interior:

\tau=(1-x)n\sigma_{0}r

\sigma_{0} the photoelectric cross section

x=n_{e}/n is the ionization fraction

The thermal energy of the gas per unit mass

E=(3/2)\kappa T(n+n_{e})/\rho

\frac{\partial E}{\partial t}=-(\bold{U}_{rel}\cdot\nabla)E-\frac{P}{\rho}(\nabla\cdot \bold{u})-\frac{\Lambda-\Gamma}{\rho}

1/n\sigma is the UV photon mean free path \tau equals the number of mean free paths in the distance r.

Data has computational units of the density of neutral hydrogen, ionized hydrogen, ionization rate, and heating rate
tau =sigmaH*nH+sigmad*nHII

Discretization

Discretize the equation of n_{e} in the same way as ​Thermal Conduction

\partial_{t}n_{e} +\partial_{i}(n_{e}u_{i})=I-r

Replace n_{e} with n_{e}^{*} where
n_{e}^{*}=(1-\Phi)n_{e}+\Phi n_{e}^{'}
and \Phi=1 for backward Euler, \Phi=1/2 for Crank Nicholson

\partial_{t}n_{e} +u_{i}\partial_{i}n_{e}^{*}+n_{e}^{*}\partial_{i}u_{i}=I-r

Links and References

1. ​free-streaming limit
   

2. ​Flux limiter

3. ​Thermal Conduction

4. ​Stromgren sphere problem

5. ​Ionization front physics

6. ​Yorke 1986

7. ​Radiation pressure

8. ​Radiation pressure 2

9. ​Ionization front

10. ​ionizaed

related code

tau	optical depth tau=sigmaH*nH + sigmad*nHII
scale
emtau	exp(-tau)
value range
ionization rate	flux*(1-exp(-tau))/dx	has computationa units of number per volume per second

f (lPencilLayout) then
       CALL UnloadFieldFromLayout(LineTransferLayout, data, LineTransferFields, (/.false.,.false.,.false./), 0)
    else
       Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),LineTransferFields)=data
    end if

! 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

Test

1. LineTransferTest Module

​global.data;​physics.data;​problem.data;

IonizingFlux=2d2	​movie
IonizingFlux=2d12	​movie