mod_afld.t

Go to the documentation of this file.

!> Module for including anisotropic flux limited diffusion (AFLD)-approximation in Radiation-hydrodynamics simulations using mod_rhd
!> Based on Turner and stone 2001. See
!> [1]Moens, N., Sundqvist, J. O., El Mellah, I., Poniatowski, L., Teunissen, J., and Keppens, R.,
!> Radiation-hydrodynamics with MPI-AMRVAC . Flux-limited diffusion
!> <i>Astronomy and Astrophysics</i>, vol. 657, 2022. doi:10.1051/0004-6361/202141023.
!> For more information.
 
module mod_afld
    use mod_comm_lib, only: mpistop
    implicit none
 
    !> source split for energy interact and radforce:
    logical :: fld_eint_split = .false.
    logical :: fld_radforce_split = .false.
    !> Opacity per unit of unit_density
    double precision, public :: fld_kappa0 = 0.34d0
    !> mean particle mass
    double precision, public :: afld_mu = 0.6d0
    !> Tolerance for bisection method for Energy sourceterms
    !> This is a percentage of the minimum of gas- and radiation energy
    double precision, public :: fld_bisect_tol = 1.d-4
    !> Tolerance for adi method for radiative Energy diffusion
    double precision, public :: fld_diff_tol = 1.d-4
    !> Number for splitting the diffusion module
    double precision, public :: diff_crit
    !> Use constant Opacity?
    character(len=8), allocatable :: fld_opacity_law(:)
    character(len=50) :: fld_opal_table = 'Y09800'
    !> Diffusion limit lambda = 0.33
    character(len=16) :: fld_fluxlimiter = 'Pomraning'
    ! !> diffusion coefficient for multigrid method
    integer, allocatable :: i_diff_mg(:)
    !> Which method to find the root for the energy interaction polynomial
    character(len=8) :: fld_interaction_method = 'Halley'
    !> Set Diffusion coefficient to unity
    logical :: fld_diff_testcase = .false.
    !> Take running average for Diffusion coefficient
    logical :: diff_coef_filter = .false.
    integer :: size_d_filter = 1
    !> Take a running average over the fluxlimiter
    logical :: flux_lim_filter = .false.
    integer :: size_l_filter = 1
    !> Use or dont use lineforce opacities
    logical :: lineforce_opacities = .false.
    !> Resume run when multigrid returns error
    logical :: diffcrash_resume = .true.
    !> A copy of rhd_Gamma
    double precision, private, protected :: afld_gamma
    !> running timestep for diffusion solver, initialised as zero
    double precision :: dt_diff = 0.d0
    !> public methods
    !> these are called in mod_rhd_phys
    public :: get_afld_rad_force
    public :: get_afld_energy_interact
    public :: afld_radforce_get_dt
    public :: afld_init
    public :: afld_get_radflux
    public :: afld_get_radpress
    public :: afld_get_fluxlimiter
    public :: afld_get_opacity
  contains
 
  !> Reading in fld-list parameters from .par file
  subroutine afld_params_read(files)
    use mod_global_parameters, only: unitpar
    use mod_constants
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    namelist /fld_list/ fld_kappa0, fld_eint_split, fld_radforce_split, &
    fld_bisect_tol, fld_diff_testcase, fld_diff_tol, fld_opal_table, &
    fld_opacity_law, fld_fluxlimiter, fld_interaction_method, &
    diff_coef_filter, size_d_filter, flux_lim_filter, size_l_filter, &
    lineforce_opacities, diffcrash_resume
 
    do n = 1, size(files)
       open(unitpar, file=trim(files(n)), status="old")
       read(unitpar, fld_list, end=111)
       111 close(unitpar)
    end do
  end subroutine afld_params_read
 
  !> Initialising FLD-module:
  !> Read opacities
  !> Initialise Multigrid
  !> adimensionalise kappa
  !> Add extra variables to w-array, flux, kappa, eddington Tensor
  !> Lambda and R
  !> ...
  subroutine afld_init(He_abundance, afld_gamma)
    use mod_global_parameters
    use mod_variables
    use mod_physics
    use mod_opal_opacity, only: init_opal_table
    use mod_multigrid_coupling
    double precision, intent(in) :: he_abundance, afld_gamma
    double precision :: sigma_thomson
    integer :: idir,jdir
    character(len=1) :: ind_1
    character(len=1) :: ind_2
    character(len=2) :: cmp_f
    character(len=5) :: cmp_e
 
    !> allocate boundary conditions
    allocate(fld_opacity_law(1:ndim))
    fld_opacity_law(1:ndim) = 'const'
 
    !> read par files
    call afld_params_read(par_files)
 
    use_multigrid = .true.
    phys_implicit_update => diffuse_e_rad_mg
    phys_evaluate_implicit => evaluate_e_rad_mg
    mg%n_extra_vars = 1
    mg%operator_type = mg_ahelmholtz
    allocate(i_diff_mg(ndim))
    do idir = 1,ndim
      write(ind_1,'(I1)') idir
      cmp_f = 'D' // ind_1
      i_diff_mg(idir) = var_set_extravar(cmp_f,cmp_f)
    enddo
    !> Need mean molecular weight
    afld_mu = (1.+4*he_abundance)/(2.+3.*he_abundance)
  end subroutine afld_init
 
  !> w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
  !> This subroutine handles the radiation force
  subroutine get_afld_rad_force(qdt,ixI^L,ixO^L,wCT,w,x,&
       energy,qsourcesplit,active)
    use mod_constants
    use mod_global_parameters
    use mod_usr_methods
    use mod_geometry
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: qdt, x(ixi^s,1:ndim)
    double precision, intent(in)    :: wct(ixi^s,1:nw)
    double precision, intent(inout) :: w(ixi^s,1:nw)
    logical, intent(in) :: energy,qsourcesplit
    logical, intent(inout) :: active
    double precision :: radiation_forcect(ixo^s,1:ndim)
    double precision :: kappact(ixo^s,1:ndim)
    double precision :: rad_fluxct(ixo^s,1:ndim)
    double precision :: div_v(ixi^s,1:ndim,1:ndim)
    double precision :: edd(ixo^s,1:ndim,1:ndim)
    double precision :: nabla_vp(ixo^s)
    double precision :: vel(ixi^s), grad_v(ixi^s), grad0_v(ixo^s)
    double precision :: grad_e(ixi^s)
    integer :: idir, jdir
    double precision :: afld_r(ixo^s,1:ndim), lambda(ixo^s,1:ndim)
 
    !> Calculate and add sourceterms
    if(qsourcesplit .eqv. fld_radforce_split) then
      active = .true.
      call afld_get_fluxlimiter(wct,x,ixi^l,ixo^l,lambda,afld_r,nghostcells)
      do idir = 1,ndim
        !> Radiation force = kappa*rho/c *Flux = lambda gradE
        call gradient(wct(ixi^s,iw_r_e),ixi^l,ixo^l,idir,grad_e,nghostcells)
        radiation_forcect(ixo^s,idir) = -lambda(ixo^s,idir)*grad_e(ixo^s)
        !> Momentum equation source term
        w(ixo^s,iw_mom(idir)) = w(ixo^s,iw_mom(idir)) &
            + qdt*radiation_forcect(ixo^s,idir)
        !> Energy equation source term (kinetic energy)
        w(ixo^s,iw_e) = w(ixo^s,iw_e) &
            + qdt*wct(ixo^s,iw_mom(idir))/wct(ixo^s,iw_rho)*radiation_forcect(ixo^s,idir)
      enddo
      !> Photon tiring
      !> calculate tensor div_v
      !> !$OMP PARALLEL DO
      do idir = 1,ndim
        do jdir = 1,ndim
          vel(ixi^s) = wct(ixi^s,iw_mom(jdir))/wct(ixi^s,iw_rho)
          call gradient(vel,ixi^l,ixo^l,idir,grad_v)
          div_v(ixo^s,idir,jdir) = grad_v(ixo^s)
        enddo
      enddo
      !> !$OMP END PARALLEL DO
      call afld_get_eddington(wct, x, ixi^l, ixo^l, edd, nghostcells)
      !> VARIABLE NAMES DIV ARE ACTUALLY GRADIENTS
      {^ifoned
      nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1) 
      }
      {^iftwod
      !>eq 34 Turner and stone (Only 2D)
      nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1)  &
                      + div_v(ixo^s,1,2)*edd(ixo^s,1,2)  &
                      + div_v(ixo^s,2,1)*edd(ixo^s,2,1)  &
                      + div_v(ixo^s,2,2)*edd(ixo^s,2,2)
      }
      {^ifthreed
      nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1)  &
                      + div_v(ixo^s,1,2)*edd(ixo^s,1,2)  &
                      + div_v(ixo^s,1,3)*edd(ixo^s,1,3)  &
                      + div_v(ixo^s,2,1)*edd(ixo^s,2,1)  &
                      + div_v(ixo^s,2,2)*edd(ixo^s,2,2)  &
                      + div_v(ixo^s,2,3)*edd(ixo^s,2,3)  &
                      + div_v(ixo^s,3,1)*edd(ixo^s,3,1)  &
                      + div_v(ixo^s,3,2)*edd(ixo^s,3,2)  &
                      + div_v(ixo^s,3,3)*edd(ixo^s,3,3)
      }
      w(ixo^s,iw_r_e) = w(ixo^s,iw_r_e) &
          - qdt * nabla_vp(ixo^s)*wct(ixo^s,iw_r_e)
    end if
  end subroutine get_afld_rad_force
 
  subroutine afld_radforce_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
    use mod_global_parameters
    use mod_usr_methods
    use mod_geometry
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: dx^d, x(ixi^s,1:ndim), w(ixi^s,1:nw)
    double precision, intent(inout) :: dtnew
    integer :: idir
    double precision :: radiation_force(ixo^s,1:ndim)
    double precision :: lambda(ixo^s,1:ndim), grad_e(ixi^s)
    double precision :: dxinv(1:ndim), afld_r(ixo^s,1:ndim)
    double precision :: max_grad
 
    ^d&dxinv(^d)=one/dx^d;
    call afld_get_fluxlimiter(w,x,ixi^l,ixo^l,lambda,afld_r,nghostcells)
    do idir = 1, ndim
      call gradient(w(ixi^s,iw_r_e),ixi^l,ixo^l,idir,grad_e,nghostcells)
      radiation_force(ixo^s,idir) = -lambda(ixo^s,idir)*grad_e(ixo^s)
      max_grad = maxval(abs(radiation_force(ixo^s,idir)))
      max_grad = max(max_grad, epsilon(1.0d0))
      dtnew = min(dtnew, courantpar/sqrt(max_grad*dxinv(idir)))
    end do
  end subroutine afld_radforce_get_dt
 
  !> w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
  !> This subroutine handles the energy exchange between gas and radiation
  subroutine get_afld_energy_interact(qdt,ixI^L,ixO^L,wCT,w,x,&
       energy,qsourcesplit,active)
    use mod_constants
    use mod_global_parameters
    use mod_usr_methods
    use mod_physics, only: phys_get_pthermal !needed to get temp
    integer, intent(in)             :: ixi^l, ixo^l
    double precision, intent(in)    :: qdt, x(ixi^s,1:ndim)
    double precision, intent(in)    :: wct(ixi^s,1:nw)
    double precision, intent(inout) :: w(ixi^s,1:nw)
    logical, intent(in) :: energy,qsourcesplit
    logical, intent(inout) :: active
 
    !> Calculate and add sourceterms
    if(qsourcesplit .eqv. fld_eint_split) then
      active = .true.
      !> Add energy sourceterms
      call energy_interaction(w,w,x,ixi^l,ixo^l)
    end if
  end subroutine get_afld_energy_interact
 
  !> Sets the opacity in the w-array
  !> by calling mod_opal_opacity
  subroutine afld_get_opacity(w, x, ixI^L, ixO^L, fld_kappa)
    use mod_global_parameters
    use mod_physics, only: phys_get_pthermal
    use mod_physics, only: phys_get_tgas
    use mod_usr_methods
    use mod_opal_opacity
    integer, intent(in)           :: ixi^l, ixo^l
    double precision, intent(in)  :: w(ixi^s, 1:nw)
    double precision, intent(in)  :: x(ixi^s, 1:ndim)
    double precision, intent(out) :: fld_kappa(ixo^s,1:ndim)
    integer :: i,j,ix^d, idir
    double precision :: temp(ixi^s), pth(ixi^s), a2(ixo^s)
    double precision :: rho0,temp0,n,sigma_b
    double precision :: akram, bkram
    double precision :: vth(ixo^s), gradv(ixi^s), eta(ixo^s), t(ixo^s)
 
    do idir = 1, ndim
      select case (fld_opacity_law(idir))
      case('const')
        fld_kappa(ixo^s,idir) = fld_kappa0/unit_opacity
      case('thomson')
        fld_kappa(ixo^s,idir) = fld_kappa0/unit_opacity
      case('kramers')
        rho0 = half !> Take lower value of rho in domain
        fld_kappa(ixo^s,idir) = fld_kappa0/unit_opacity*((w(ixo^s,iw_rho)/rho0))
      case('bump')
        !> Opacity bump
        rho0 = 0.2d0
        n = 7.d0
        sigma_b = 2.d-2
        fld_kappa(ixo^s,idir) = fld_kappa0/unit_opacity*(one + n*dexp(-one/sigma_b*(dlog(w(ixo^s,iw_rho)/rho0))**two))
      case('non_iso')
        call phys_get_pthermal(w,x,ixi^l,ixo^l,temp)
        temp(ixo^s) = temp(ixo^s)/w(ixo^s,iw_rho)
        rho0 = 0.5d0 !> Take lower value of rho in domain
        temp0 = one
        n = -7.d0/two
        fld_kappa(ixo^s,idir) = fld_kappa0/unit_opacity*(w(ixo^s,iw_rho)/rho0)*(temp(ixo^s)/temp0)**n
      case('fastwind')
        call phys_get_pthermal(w,x,ixi^l,ixo^l,pth)
        a2(ixo^s) = pth(ixo^s)/w(ixo^s,iw_rho)*unit_velocity**2.d0
        akram = 13.1351597305d0
        bkram = -4.5182188206d0
        fld_kappa(ixo^s,idir) = fld_kappa0/unit_opacity &
            * (1.d0+10.d0**akram*w(ixo^s,iw_rho)*unit_density*(a2(ixo^s)/1.d12)**bkram)
        {do ix^d=ixomin^d,ixomax^d\ }
          !> Hard limit on kappa
          fld_kappa(ix^d,idir) = min(fld_kappa(ix^d,idir),2.3d0*fld_kappa0/unit_opacity)
        {enddo\}
      case('opal')
        call phys_get_tgas(w,x,ixi^l,ixo^l,temp)
        {do ix^d=ixomin^d,ixomax^d\ }
          rho0 = w(ix^d,iw_rho)*unit_density
          temp0 = temp(ix^d)*unit_temperature
          call set_opal_opacity(rho0,temp0,n)
          fld_kappa(ix^d,idir) = n/unit_opacity
        {enddo\}
      case('special')
        if (.not. associated(usr_special_aniso_opacity)) then
          call mpistop("special opacity not defined")
        endif
        call usr_special_aniso_opacity(ixi^l,ixo^l,w,x,fld_kappa(ixo^s,idir),idir)
      case default
        call mpistop("Doesn't know opacity law")
      end select
    end do
  end subroutine afld_get_opacity
 
  !> Calculate fld flux limiter
  !> This subroutine calculates flux limiter lambda using the prescription
  !> stored in fld_fluxlimiter.
  !> It also calculates the ratio of radiation scaleheight and mean free path
  subroutine afld_get_fluxlimiter(w,x,ixI^L,ixO^L,fld_lambda,fld_R,nth)
    use mod_global_parameters
    use mod_geometry
    use mod_usr_methods
    integer, intent(in)          :: ixi^l, ixo^l, nth
    double precision, intent(in) :: w(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    double precision, intent(out) :: fld_r(ixo^s,1:ndim), fld_lambda(ixo^s,1:ndim)
    integer :: idir, ix^d
    integer :: filter, idim
    double precision :: kappa(ixo^s,1:ndim)
    double precision ::  normgrad2(ixo^s)
    double precision :: grad_r_e(ixi^s)
    double precision :: rad_e(ixi^s)
    double precision :: tmp_l(ixi^s), filtered_l(ixi^s)
 
    select case (fld_fluxlimiter)
    case('Diffusion')
      fld_lambda(ixo^s,1:ndim) = 1.d0/3.d0
      fld_r(ixo^s,1:ndim) = zero
    case('FreeStream')
      !> Calculate R everywhere
      !> |grad E|/(rho kappa E)
      rad_e(ixi^s) = w(ixi^s, iw_r_e)
      call afld_get_opacity(w,x,ixi^l,ixo^l,kappa)
      do idir=1,ndim
        call gradient(rad_e,ixi^l,ixo^l,idir,grad_r_e,nth)
        normgrad2(ixo^s) = grad_r_e(ixo^s)**2
        fld_r(ixo^s,idir) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s,idir)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
        ! Note: possible zeros for uniform setup, to adjust
        if(maxval(normgrad2(ixo^s))==0.0d0)call mpistop("choose other fluxlimiter, FreeStream has issues")
        !> Calculate the flux limiter, lambda
        fld_lambda(ixo^s,idir) = one/fld_r(ixo^s,idir)
      end do
    case('Pomraning')
      !> Calculate R everywhere
      !> |grad E|/(rho kappa E)
      normgrad2(ixo^s) = zero
      rad_e(ixi^s) = w(ixi^s,iw_r_e)
      call afld_get_opacity(w,x,ixi^l,ixo^l,kappa)
      do idir = 1,ndim
        call gradient(rad_e,ixi^l,ixo^l,idir,grad_r_e,nth)
        normgrad2(ixo^s) = grad_r_e(ixo^s)**2
        fld_r(ixo^s,idir) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s,idir)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
        !> Calculate the flux limiter, lambda
        !> Levermore and Pomraning: lambda = (2 + R)/(6 + 3R + R^2)
        fld_lambda(ixo^s,idir) = (2.d0+fld_r(ixo^s,idir))/(6.d0+3*fld_r(ixo^s,idir)+fld_r(ixo^s,idir)**2.d0)
      end do
    case('Minerbo')
      !> Calculate R everywhere
      !> |grad E|/(rho kappa E)
      normgrad2(ixo^s) = zero
      rad_e(ixi^s) = w(ixi^s,iw_r_e)
      call afld_get_opacity(w,x,ixi^l,ixo^l,kappa)
      do idir = 1,ndim
        call gradient(rad_e,ixi^l,ixo^l,idir,grad_r_e,nth)
        normgrad2(ixo^s) = grad_r_e(ixo^s)**2
        fld_r(ixo^s,idir) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s,idir)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
        !> Calculate the flux limiter, lambda
        !> Minerbo:
        {do ix^d=ixomin^d,ixomax^d\ }
        if(fld_r(ix^d,idir) .lt. 3.d0/2.d0) then
          fld_lambda(ix^d,idir) = 2.d0/(3.d0 + dsqrt(9.d0 + 12.d0*fld_r(ix^d,idir)**2.d0))
        else
          fld_lambda(ix^d,idir) = 1.d0/(1.d0 + fld_r(ix^d,idir) + dsqrt(1.d0 + 2.d0*fld_r(ix^d,idir)))
        endif
        {enddo\}
      end do
    case('special')
      if (.not. associated(usr_special_fluxlimiter)) then
        call mpistop("special fluxlimiter not defined")
      endif
      call usr_special_fluxlimiter(ixi^l, ixo^l, w, x, fld_lambda, fld_r)
    case default
      call mpistop('Fluxlimiter unknown')
    end select
 
    if(flux_lim_filter) then
      if(size_l_filter .lt. 1) call mpistop("D filter of size < 1 makes no sense")
      if(size_l_filter .gt. nghostcells) call mpistop("D filter of size > nghostcells makes no sense")
      do idir=1,ndim
        tmp_l(ixo^s) = fld_lambda(ixo^s,idir)
        filtered_l(ixo^s) = zero
        do filter = 1,size_l_filter
          {do ix^d = ixomin^d+size_d_filter,ixomax^d-size_l_filter\}
            do idim = 1,ndim
              filtered_l(ix^d) = filtered_l(ix^d) &
                               + tmp_l(ix^d+filter*kr(idim,^d)) &
                               + tmp_l(ix^d-filter*kr(idim,^d))
            enddo
          {enddo\}
        enddo
        {do ix^d = ixomin^d+size_d_filter,ixomax^d-size_d_filter\}
          tmp_l(ix^d) = (tmp_l(ix^d)+filtered_l(ix^d))/(1+2*size_l_filter*ndim)
        {enddo\}
      enddo
      fld_lambda(ixo^s,idir) = tmp_l(ixo^s)
    endif
  end subroutine afld_get_fluxlimiter
 
  !> Calculate Radiation Flux
  !> Calculate the Flux using the fld closure relation
  !> F = -c*lambda/(kappa*rho) *grad E
  subroutine afld_get_radflux(w, x, ixI^L, ixO^L, rad_flux, nth)
    use mod_global_parameters
    use mod_geometry
    integer, intent(in)          :: ixi^l, ixo^l, nth
    double precision, intent(in) :: w(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    double precision, intent(out) :: rad_flux(ixo^s, 1:ndim)
    integer :: ix^d, idir
    double precision :: cc
    double precision :: grad_r_e(ixi^s)
    double precision :: rad_e(ixi^s)
    double precision :: kappa(ixo^s,1:ndim), lambda(ixo^s,1:ndim), fld_r(ixo^s,1:ndim)
 
    cc = const_c/unit_velocity
    rad_e(ixi^s) = w(ixi^s, iw_r_e)
    call afld_get_opacity(w, x, ixi^l, ixo^l, kappa)
    call afld_get_fluxlimiter(w, x, ixi^l, ixo^l, lambda, fld_r, nghostcells)
    do idir = 1,ndim
      call gradient(rad_e,ixi^l,ixo^l,idir,grad_r_e,nghostcells)
      rad_flux(ixo^s, idir) = -cc*lambda(ixo^s,idir)/(kappa(ixo^s,idir)*w(ixo^s,iw_rho))*grad_r_e(ixo^s)
    end do
  end subroutine afld_get_radflux
 
  !> Calculate Eddington-tensor
  !> Stores Eddington-tensor in w-array
  subroutine afld_get_eddington(w, x, ixI^L, ixO^L, eddington_tensor, nth)
    use mod_global_parameters
    use mod_geometry
    integer, intent(in)          :: ixI^L, ixO^L, nth
    double precision, intent(in) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision, intent(out) :: eddington_tensor(ixO^S,1:ndim,1:ndim)
    double precision :: tnsr2(ixO^S,1:ndim,1:ndim)
    double precision :: normgrad2(ixO^S), f(ixO^S,1:ndim)
    double precision :: grad_r_e(ixI^S,1:ndim), rad_e(ixI^S)
    double precision :: lambda(ixO^S,1:ndim), fld_R(ixO^S,1:ndim)
    integer :: i,j, idir,jdir
 
    call afld_get_fluxlimiter(w,x,ixi^l,ixo^l,lambda,fld_r, nth)
    !> Calculate R everywhere
    !> |grad E|/(rho kappa E)
    normgrad2(ixo^s) = zero
    rad_e(ixi^s) = w(ixi^s, iw_r_e)
    do idir = 1,ndim
      call gradient(rad_e,ixi^l,ixo^l,idir,grad_r_e(ixi^s,idir),nth)
      normgrad2(ixo^s) = normgrad2(ixo^s)+grad_r_e(ixo^s,idir)**2
      !> Calculate radiation pressure
      !> P = (lambda + lambda^2 R^2)*E
      f(ixo^s,idir) = lambda(ixo^s,idir)+lambda(ixo^s,idir)**two*fld_r(ixo^s,idir)**two
      f(ixo^s,idir) = half*(one-f(ixo^s,idir))+half*(3.d0*f(ixo^s,idir)-one)
    end do
 
    {^ifoned
    eddington_tensor(ixo^s,1,1) = f(ixo^s,1)
    }
    {^nooned
    do idir = 1,ndim
      eddington_tensor(ixo^s,idir,idir) = half*(one-f(ixo^s,idir))
    enddo
    do idir = 1,ndim
      do jdir = 1,ndim
        if(idir .ne. jdir) eddington_tensor(ixo^s,idir,jdir) = zero
        tnsr2(ixo^s,idir,jdir) =  half*(3.d0*f(ixo^s,idir) - 1)&
          * grad_r_e(ixo^s,idir)*grad_r_e(ixo^s,jdir)/max(normgrad2(ixo^s),smalldouble)
      enddo
    enddo
    do idir = 1,ndim
      do jdir = 1,ndim
        where(normgrad2(ixo^s) .gt. smalldouble)
          eddington_tensor(ixo^s,idir,jdir) = eddington_tensor(ixo^s,idir,jdir)+tnsr2(ixo^s,idir,jdir)
        endwhere
      enddo
    enddo
    }
  end subroutine afld_get_eddington
 
  !> Calculate Radiation Pressure
  !> Returns Radiation Pressure as tensor
  subroutine afld_get_radpress(w, x, ixI^L, ixO^L, rad_pressure, nth)
    use mod_global_parameters
    integer, intent(in)          :: ixi^l, ixo^l, nth
    double precision, intent(in) :: w(ixi^s, 1:nw)
    double precision, intent(in) :: x(ixi^s, 1:ndim)
    double precision, intent(out):: rad_pressure(ixo^s,1:ndim,1:ndim)
    integer :: i,j
    double precision :: eddington_tensor(ixo^s,1:ndim,1:ndim)
 
    call afld_get_eddington(w,x,ixi^l,ixo^l,eddington_tensor, nth)
    do i=1,ndim
      do j=1,ndim
        rad_pressure(ixo^s,i,j) = eddington_tensor(ixo^s,i,j)*w(ixo^s,iw_r_e)
      enddo
    enddo
  end subroutine afld_get_radpress
 
!> Calling all subroutines to perform the multigrid method
!> Communicates rad_e and diff_coeff to multigrid library
  subroutine diffuse_e_rad_mg(dtfactor,qdt,qtC,psa,psb)
    use mod_global_parameters
    use mod_forest
    use mod_ghostcells_update
    use mod_multigrid_coupling
    use mod_physics, only: phys_set_mg_bounds
    type(state), target :: psa(max_blocks)
    type(state), target :: psb(max_blocks)
    double precision, intent(in) :: qdt
    double precision, intent(in) :: qtC
    double precision, intent(in) :: dtfactor
    integer, parameter           :: max_its = 50
    integer                      :: n
    integer                      :: iw_to,iw_from
    integer                      :: iigrid, igrid, id
    integer                      :: nc, lvl, i
    type(tree_node), pointer     :: pnode
    double precision             :: fac, facD
    double precision             :: res, max_residual, lambda
 
    !> Set diffusion timestep, add previous timestep if mg did not converge:
    dt_diff = dt_diff + qdt
    ! Avoid setting a very restrictive limit to the residual when the time step
    ! is small (as the operator is ~ 1/(D * qdt))
    if (qdt < dtmin) then
        if(mype==0)then
            print *,'skipping implicit solve: dt too small!'
            print *,'Currently at time=',global_time,' time step=',qdt,' dtmin=',dtmin
        endif
        return
    endif
    ! max_residual = 1d-7/qdt
    max_residual = fld_diff_tol !1d-7/qdt
    mg%operator_type = mg_ahelmholtz
    mg%smoother_type = mg_smoother_gsrb
    call mg_set_methods(mg)
    if (.not. mg%is_allocated) call mpistop("multigrid tree not allocated yet")
    lambda = 1.d0/(dtfactor * dt_diff)
    call ahelmholtz_set_lambda(lambda)
    call update_diffcoeff(psb)
    fac = 1.d0
    facd = 1.d0
    !This is mg_copy_to_tree from psb state
    {^ifoned
    call mg_copy_to_tree(i_diff_mg(1), mg_iveps, factor=facd, state_from=psb)
    if (time_advance)then
      call mg_restrict(mg, mg_iveps)
      call mg_fill_ghost_cells(mg, mg_iveps)
    endif
    }
    {^iftwod
    call mg_copy_to_tree(i_diff_mg(1), mg_iveps1, factor=facd, state_from=psb)
    call mg_copy_to_tree(i_diff_mg(2), mg_iveps2, factor=facd, state_from=psb)
    if (time_advance)then
      call mg_restrict(mg, mg_iveps1)
      call mg_fill_ghost_cells(mg, mg_iveps1)
 
      call mg_restrict(mg, mg_iveps2)
      call mg_fill_ghost_cells(mg, mg_iveps2)
    endif
    }
    !This is mg_copy_to_tree from psb state
    call mg_copy_to_tree(iw_r_e, mg_iphi, factor=fac, state_from=psb)
    !>replace call set_rhs(mg, -1/dt, 0.0_dp)
    call mg_copy_to_tree(iw_r_e, mg_irhs, factor=-1/(dtfactor*dt_diff), state_from=psb)
    call phys_set_mg_bounds()
    call mg_fas_fmg(mg, .true., max_res=res)
    do n = 1, max_its
      ! print*, n, res
      if (res < max_residual) exit
       call mg_fas_vcycle(mg, max_res=res)
    end do
    if (res .le. 0.d0) then
      if (diffcrash_resume) then
        if (mg%my_rank == 0) &
        write(*,*) it, ' residual zero ', res
        return
      endif
      if (mg%my_rank == 0) then
        print*, res
        error stop "Diffusion residual to zero"
      endif
    endif
    if(n == max_its + 1) then
      if(diffcrash_resume) then
        if(mg%my_rank == 0) &
        write(*,*) it, ' residual high ', res
        return
      endif
      if(mg%my_rank == 0) then
        print *, "Did you specify boundary conditions correctly?"
        print *, "Or is the variation in diffusion too large?"
        print*, n, res
        print *, mg%bc(1, mg_iphi)%bc_value, mg%bc(2, mg_iphi)%bc_value
      end if
      error stop "diffusion_solve: no convergence"
    end if
    !> Reset dt_diff when diffusion worked out
    dt_diff = 0.d0
    ! !This is mg_copy_from_tree_gc for psa state
    call mg_copy_from_tree_gc(mg_iphi, iw_r_e, state_to=psa)
  end subroutine diffuse_e_rad_mg
 
  !> inplace update of psa==>F_im(psa)
  subroutine evaluate_e_rad_mg(qtC,psa)
    use mod_global_parameters
    use mod_ghostcells_update
    type(state), target :: psa(max_blocks)
    double precision, intent(in) :: qtC
    integer :: iigrid, igrid, level
    integer :: ixO^L
 
    call update_diffcoeff(psa)
    ixo^l=ixg^ll^lsub1;
    !$OMP PARALLEL DO PRIVATE(igrid)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
       ! call afld_get_diffcoef_central(psa(igrid)%w, psa(igrid)%w, psa(igrid)%x, ixG^LL, ixO^L)
       call put_diffterm_onegrid(ixg^ll,ixo^l,psa(igrid)%w)
    end do
    !$OMP END PARALLEL DO
    ! enforce boundary conditions for psa
    call getbc(qtc,0.d0,psa,1,nwflux+nwaux)
  end subroutine evaluate_e_rad_mg
 
  !> inplace update of psa==>F_im(psa)
  subroutine put_diffterm_onegrid(ixI^L,ixO^L,w)
    use mod_global_parameters
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(inout) :: w(ixI^S, 1:nw)
    double precision :: gradE(ixO^S),divF(ixO^S)
    double precision :: divF_h(ixO^S),divF_j(ixO^S)
    double precision :: diff_term(ixO^S)
    integer                       :: idir, jxO^L, hxO^L
 
    divf(ixo^s) = 0.d0
    do idir = 1,ndim
      hxo^l=ixo^l-kr(idir,^d);
      jxo^l=ixo^l+kr(idir,^d);
      divf_h(ixo^s) = w(ixo^s,i_diff_mg(idir))*w(hxo^s,i_diff_mg(idir))/(w(ixo^s,i_diff_mg(idir)) + w(hxo^s,i_diff_mg(idir)))*(w(hxo^s,iw_r_e) - w(ixo^s,iw_r_e))
      divf_j(ixo^s) = w(ixo^s,i_diff_mg(idir))*w(jxo^s,i_diff_mg(idir))/(w(ixo^s,i_diff_mg(idir)) + w(jxo^s,i_diff_mg(idir)))*(w(jxo^s,iw_r_e) - w(ixo^s,iw_r_e))
      divf(ixo^s) = divf(ixo^s) + 2.d0*(divf_h(ixo^s) + divf_j(ixo^s))/dxlevel(idir)**2
    enddo
    w(ixo^s,iw_r_e) = divf(ixo^s)
  end subroutine put_diffterm_onegrid
 
  !> Calculates cell-centered diffusion coefficient to be used in multigrid
  subroutine afld_get_diffcoef_central(w, wCT, x, ixI^L, ixO^L)
    use mod_global_parameters
    use mod_geometry
    use mod_usr_methods
    integer, intent(in)          :: ixI^L, ixO^L
    double precision, intent(inout) :: w(ixI^S, 1:nw)
    double precision, intent(in) :: wCT(ixI^S, 1:nw)
    double precision, intent(in) :: x(ixI^S, 1:ndim)
    double precision :: kappa(ixO^S,1:ndim), lambda(ixO^S,1:ndim), fld_R(ixO^S,1:ndim)
    double precision :: max_D(ixI^S), grad_r_e(ixI^S), rad_e(ixI^S)
    double precision :: cc
    integer :: idir,i,j,ix^D
 
    cc = const_c/unit_velocity
    if(fld_diff_testcase) then
      w(ixo^s,i_diff_mg(:)) = 1.d0
    else
      call afld_get_opacity(wct, x, ixi^l, ixo^l, kappa)
      call afld_get_fluxlimiter(wct, x, ixi^l, ixo^l, lambda, fld_r, nghostcells-1)
      do idir = 1,ndim
        !> calculate diffusion coefficient
        w(ixo^s,i_diff_mg(idir)) = cc*lambda(ixo^s,idir)/(kappa(ixo^s,idir)*wct(ixo^s,iw_rho))
        where(w(ixo^s,i_diff_mg(idir)) .lt. 0.d0)
          w(ixo^s,i_diff_mg(idir)) = smalldouble
        end where
        if(diff_coef_filter) then
          call afld_smooth_diffcoef(w,ixi^l,ixo^l,idir)
        endif
      enddo
    endif
    if(associated(usr_special_diffcoef)) &
      call usr_special_diffcoef(w, wct, x, ixi^l, ixo^l)
  end subroutine afld_get_diffcoef_central
 
  subroutine update_diffcoeff(psa)
    use mod_global_parameters
    type(state), target :: psa(max_blocks)
    integer :: iigrid, igrid, level
    integer :: ixO^L
 
    ixo^l=ixg^ll^lsub1;
    !$OMP PARALLEL DO PRIVATE(igrid)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
        call afld_get_diffcoef_central(psa(igrid)%w, psa(igrid)%w, psa(igrid)%x, ixg^ll, ixo^l)
    end do
    !$OMP END PARALLEL DO
  end subroutine update_diffcoeff
 
  !> Use running average on Diffusion coefficient
  subroutine afld_smooth_diffcoef(w, ixI^L, ixO^L,idir)
    use mod_global_parameters
    integer, intent(in) :: ixI^L, ixO^L, idir
    double precision, intent(inout) :: w(ixI^S, 1:nw)
    double precision :: tmp_D(ixI^S), filtered_D(ixI^S)
    integer :: ix^D, filter, idim
 
    if(size_d_filter .lt. 1) call mpistop("D filter of size < 1 makes no sense")
    if(size_d_filter .gt. nghostcells) call mpistop("D filter of size > nghostcells makes no sense")
    tmp_d(ixo^s) = w(ixo^s,i_diff_mg(idir))
    filtered_d(ixo^s) = zero
    do filter = 1,size_d_filter
      {do ix^d = ixomin^d+size_d_filter,ixomax^d-size_d_filter\}
        do idim = 1,ndim
          filtered_d(ix^d) = filtered_d(ix^d) &
                           + tmp_d(ix^d+filter*kr(idim,^d)) &
                           + tmp_d(ix^d-filter*kr(idim,^d))
        enddo
      {enddo\}
    enddo
    {do ix^d = ixomin^d+size_d_filter,ixomax^d-size_d_filter\}
      tmp_d(ix^d) = (tmp_d(ix^d)+filtered_d(ix^d))/(1+2*size_d_filter*ndim)
    {enddo\}
    w(ixo^s,i_diff_mg(idir)) = tmp_d(ixo^s)
  end subroutine afld_smooth_diffcoef
 
  !> This subroutine calculates the radiation heating, radiation cooling
  !> and photon tiring using an implicit scheme.
  !> These sourceterms are applied using the root-finding of a 4th order polynomial
  !> This routine loops over every cell in the domain
  !> and computes the coefficients of the polynomials in every cell
  subroutine energy_interaction(w, wCT, x, ixI^L, ixO^L)
    use mod_global_parameters
    use mod_geometry
    use mod_physics
    use mod_usr_methods
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: x(ixI^S,1:ndim)
    double precision, intent(in)    :: wCT(ixI^S,1:nw)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    double precision :: a1(ixO^S), a2(ixO^S)
    double precision :: c0(ixO^S), c1(ixO^S)
    double precision :: e_gas(ixO^S), E_rad(ixO^S)
    double precision :: kappa(ixO^S)
    double precision :: sigma_b,cc
    integer :: i,j,ix^D
 
    cc = const_c/unit_velocity
    sigma_b = const_rad_a*cc/4.d0*(unit_temperature**4.d0)/(unit_pressure)
    !> e_gas is the INTERNAL ENERGY without KINETIC ENERGY
    e_gas(ixo^s) = wct(ixo^s,iw_e)-half*sum(wct(ixo^s, iw_mom(:))**2, dim=ndim+1)/wct(ixo^s, iw_rho)
    if(allocated(iw_mag)) then
      e_gas(ixo^s) = e_gas(ixo^s)-half*sum(wct(ixo^s,iw_mag(:))**2,dim=ndim+1)
    endif
    {do ix^d=ixomin^d,ixomax^d\ }
      e_gas(ix^d) = max(e_gas(ix^d),small_pressure/(afld_gamma-1.d0))
    {enddo\}
    e_rad(ixo^s) = wct(ixo^s,iw_r_e)
    if(associated(usr_special_opacity_qdot)) then
      call usr_special_opacity_qdot(ixi^l,ixo^l,w,x,kappa)
    else
      call afld_get_opacity(wct,x,ixi^l,ixo^l,kappa)
    endif
    if(fld_interaction_method .eq. 'Instant') then
      a1(ixo^s) = 4.d0*sigma_b/cc*(afld_gamma-one)**4.d0/wct(ixo^s,iw_rho)**4.d0
      a2(ixo^s) = e_gas(ixo^s) + e_rad(ixo^s)
      c0(ixo^s) = a2(ixo^s)/a1(ixo^s)
      c1(ixo^s) = 1.d0/a1(ixo^s)
    else
      !> Calculate coefficients for polynomial
      a1(ixo^s) = 4.d0*kappa(ixo^s)*sigma_b*(afld_gamma-one)**4.d0/wct(ixo^s,iw_rho)**3.d0*dt
      a2(ixo^s) = (const_c/unit_velocity)*kappa(ixo^s)*wct(ixo^s,iw_rho)*dt
      c0(ixo^s) = ((one+a2(ixo^s))*e_gas(ixo^s)+a2(ixo^s)*e_rad(ixo^s))/a1(ixo^s)
      c1(ixo^s) = (one+a2(ixo^s))/a1(ixo^s)
    endif
    !> Loop over every cell for rootfinding method
    {do ix^d=ixomin^d,ixomax^d\ }
      select case(fld_interaction_method)
      case('Bisect')
        call bisection_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
      case('Newton')
        call newton_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
      case('Halley')
        call halley_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
      case('Instant')
        call halley_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
      case default
        call mpistop('root-method not known')
      end select
    {enddo\}
 
    if(fld_interaction_method .eq. 'Instant') then
      e_rad(ixo^s) = a2(ixo^s)-e_gas(ixo^s)
    else
      !> advance E_rad
      e_rad(ixo^s) = (a1(ixo^s)*e_gas(ixo^s)**4.d0+e_rad(ixo^s))/(one + a2(ixo^s))
    endif
    !> Update gas-energy in w, internal + kinetic
    w(ixo^s,iw_e) = e_gas(ixo^s)
    !> Beginning of module substracted wCT Ekin
    !> So now add wCT Ekin
    w(ixo^s,iw_e) = w(ixo^s,iw_e) + half*sum(wct(ixo^s, iw_mom(:))**2,dim=ndim+1)/wct(ixo^s, iw_rho)
    if(allocated(iw_mag)) then
      w(ixo^s,iw_e) = w(ixo^s,iw_e)+half*sum(wct(ixo^s, iw_mag(:))**2,dim=ndim+1)
    endif
    !> Update rad-energy in w
    w(ixo^s,iw_r_e) = e_rad(ixo^s)
  end subroutine energy_interaction
 
  !> Find the root of the 4th degree polynomial using the bisection method
  subroutine bisection_method(e_gas, E_rad, c0, c1)
    use mod_global_parameters
    double precision, intent(in)    :: c0, c1
    double precision, intent(in)    :: E_rad
    double precision, intent(inout) :: e_gas
    double precision :: bisect_a, bisect_b, bisect_c
    integer :: n, max_its
 
    n = 0
    max_its = 10000000
    bisect_a = zero
    bisect_b = 1.2d0*max(abs(c0/c1),abs(c0)**(1.d0/4.d0))
    do while(abs(bisect_b-bisect_a) .ge. fld_bisect_tol*min(e_gas,e_rad))
      bisect_c = (bisect_a + bisect_b)/two
      n = n +1
      if(n .gt. max_its) then
        goto 2435
        call mpistop('No convergece in bisection scheme')
      endif
      if(polynomial_bisection(bisect_a, c0, c1)*&
      polynomial_bisection(bisect_b, c0, c1) .lt. zero) then
        if(polynomial_bisection(bisect_a, c0, c1)*&
        polynomial_bisection(bisect_c, c0, c1) .lt. zero) then
          bisect_b = bisect_c
        elseif(polynomial_bisection(bisect_b, c0, c1)*&
        polynomial_bisection(bisect_c, c0, c1) .lt. zero) then
          bisect_a = bisect_c
        elseif(polynomial_bisection(bisect_a, c0, c1) .eq. zero) then
          bisect_b = bisect_a
          bisect_c = bisect_a
          goto 2435
        elseif(polynomial_bisection(bisect_b, c0, c1) .eq. zero) then
          bisect_a = bisect_b
          bisect_c = bisect_b
          goto 2435
        elseif(polynomial_bisection(bisect_c, c0, c1) .eq. zero) then
          bisect_a = bisect_c
          bisect_b = bisect_c
          goto 2435
        else
          call mpistop("Problem with fld bisection method")
        endif
      elseif(polynomial_bisection(bisect_a, c0, c1) &
        - polynomial_bisection(bisect_b, c0, c1) .lt. fld_bisect_tol*polynomial_bisection(bisect_a, c0, c1)) then
        goto 2435
      else
        bisect_a = e_gas
        bisect_b = e_gas
        print*, "IGNORING GAS-RAD ENERGY EXCHANGE ", c0, c1
        print*, polynomial_bisection(bisect_a, c0, c1), polynomial_bisection(bisect_b, c0, c1)
        if(polynomial_bisection(bisect_a, c0, c1) .le. smalldouble) then
          bisect_b = bisect_a
        elseif(polynomial_bisection(bisect_a, c0, c1) .le. smalldouble) then
          bisect_a = bisect_b
        endif
        goto 2435
      endif
    enddo
      2435 e_gas = (bisect_a + bisect_b)/two
  end subroutine bisection_method
 
  !> Find the root of the 4th degree polynomial using the Newton-Ralphson method
  subroutine newton_method(e_gas, E_rad, c0, c1)
    use mod_global_parameters
    double precision, intent(in)    :: c0, c1
    double precision, intent(in)    :: E_rad
    double precision, intent(inout) :: e_gas
    double precision :: xval, yval, der, deltax
    integer :: ii
 
    yval = bigdouble
    xval = e_gas
    der = one
    deltax = one
    ii = 0
    !> Compare error with dx = dx/dy dy
    do while (abs(deltax) .gt. fld_bisect_tol)
      yval = polynomial_bisection(xval, c0, c1)
      der = dpolynomial_bisection(xval, c0, c1)
      deltax = -yval/der
      xval = xval + deltax
      ii = ii + 1
      if (ii .gt. 1d3) then
        call bisection_method(e_gas, e_rad, c0, c1)
        return
      endif
    enddo
    e_gas = xval
  end subroutine newton_method
 
  !> Find the root of the 4th degree polynomial using the Halley method
  subroutine halley_method(e_gas, E_rad, c0, c1)
    use mod_global_parameters
    double precision, intent(in)    :: c0, c1
    double precision, intent(in)    :: E_rad
    double precision, intent(inout) :: e_gas
    double precision :: xval, yval, der, dder, deltax
    integer :: ii
 
    yval = bigdouble
    xval = e_gas
    der = one
    dder = one
    deltax = one
    ii = 0
    !> Compare error with dx = dx/dy dy
    do while (abs(deltax) .gt. fld_bisect_tol)
      yval = polynomial_bisection(xval, c0, c1)
      der = dpolynomial_bisection(xval, c0, c1)
      dder = ddpolynomial_bisection(xval, c0, c1)
      deltax = -two*yval*der/(two*der**2 - yval*dder)
      xval = xval + deltax
      ii = ii + 1
      if (ii .gt. 1d3) then
        call newton_method(e_gas, e_rad, c0, c1)
        return
      endif
    enddo
    e_gas = xval
  end subroutine halley_method
 
  !> Evaluate polynomial at argument e_gas
  function polynomial_bisection(e_gas, c0, c1) result(val)
    use mod_global_parameters
    double precision, intent(in) :: e_gas
    double precision, intent(in) :: c0, c1
    double precision :: val
 
    val = e_gas**4.d0 + c1*e_gas - c0
  end function polynomial_bisection
 
  !> Evaluate first derivative of polynomial at argument e_gas
  function dpolynomial_bisection(e_gas, c0, c1) result(der)
    use mod_global_parameters
    double precision, intent(in) :: e_gas
    double precision, intent(in) :: c0, c1
    double precision :: der
 
    der = 4.d0*e_gas**3.d0 + c1
  end function dpolynomial_bisection
 
  !> Evaluate second derivative of polynomial at argument e_gas
  function ddpolynomial_bisection(e_gas, c0, c1) result(dder)
    use mod_global_parameters
    double precision, intent(in) :: e_gas
    double precision, intent(in) :: c0, c1
    double precision :: dder
 
    dder = 4.d0*3.d0*e_gas**2.d0
  end function ddpolynomial_bisection
end module mod_afld