mod_fld.t

Go to the documentation of this file.

!> Nicolas Moens with updates by RK (16/03/2026)
!> Module for including flux limited diffusion (FLD)-approximation in Radiation-hydrodynamics simulations
!> 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_fld
    use mod_comm_lib, only: mpistop
    use mod_geometry
    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 :: fld_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
    !> switches for opacity
    character(len=8)  :: fld_opacity_law = 'const'
    character(len=50) :: fld_opal_table = 'Y09800' 
    !> flux limiter choice
    character(len=16) :: fld_fluxlimiter = 'Pomraning'
    !> diffusion coefficient for multigrid method
    integer :: i_diff_mg
    !> Which method to find the root for the energy interaction polynomial
    character(len=8) :: fld_interaction_method = 'Halley'
    !> Use or dont use lineforce opacities
    logical :: lineforce_opacities = .false.
    !> Resume run when multigrid returns error
    logical :: diffcrash_resume = .true.
    !> Index for Flux weighted opacities
    integer, allocatable, public :: i_opf(:)
    !> A copy of (m)hd_Gamma
    double precision, private, protected :: fld_gamma
    !> running timestep for diffusion solver, initialised as zero
    double precision :: dt_diff = 0.d0
    !> public methods
    !> these are called in mod_hd_phys or mod_mhd_phys
    public :: fld_init
    public :: fld_get_radpress
    public :: fld_get_diffcoef_central
    public :: add_fld_rad_force
    public :: fld_radforce_get_dt
    ! these are made public for mod_usr purposes and diagnostics
    public :: fld_get_radflux
    public :: fld_get_fluxlimiter
    public :: fld_get_opacity
  contains
 
  !> Reading in fld-list parameters from .par file
  subroutine fld_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_tol,&
    fld_opacity_law, fld_fluxlimiter, fld_interaction_method, &
    lineforce_opacities, diffcrash_resume, fld_opal_table
 
    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 fld_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 fld_init(He_abundance, r_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, r_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
 
    !> read par files
    call fld_params_read(par_files)
    !> Set lineforce opacities as variable
    if(lineforce_opacities) then
      allocate(i_opf(ndim))
      do idir = 1,ndim
        write(ind_1,'(I1)') idir
        cmp_f = 'k' // ind_1
        i_opf(idir) = var_set_extravar(cmp_f,cmp_f)
      enddo
    endif
    phys_implicit_update => fld_implicit_update
    phys_evaluate_implicit => evaluate_e_rad_mg
    use_multigrid = .true.
    mg%n_extra_vars = 1
    mg%operator_type = mg_vhelmholtz
    i_diff_mg = var_set_extravar("D", "D")
    !> Need mean molecular weight
    fld_mu = (1.+4*he_abundance)/(2.+3.*he_abundance)
    !> set gamma
    fld_gamma = r_gamma
    !> Read in opacity table if necesary
    if(fld_opacity_law .eq. 'opal') call init_opal_table(fld_opal_table)
    if(fld_opacity_law .eq. 'thomson')  then
      ! special fixed value for thomson opacity, assuming cgs units
      sigma_thomson = 6.6524585d-25
      fld_kappa0 = sigma_thomson/const_mp * (1.+2.*he_abundance)/(1.+4.*he_abundance)
    endif
  end subroutine fld_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 add_fld_rad_force(qdt,ixI^L,ixO^L,wCT,wCTprim,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),wctprim(ixi^s,1:nw)
    double precision, intent(inout) :: w(ixi^s,1:nw)
    logical, intent(in) :: energy,qsourcesplit
    logical, intent(inout) :: active
    integer :: idir, jdir
    double precision :: radiation_forcect(ixo^s,1:ndim)
    double precision :: rad_fluxct(ixo^s,1:ndim)
    double precision :: grad_e(ixi^s)
    double precision :: lambda(ixo^s),fld_r(ixo^s)
 
    !> Calculate and add sourceterms
    if(qsourcesplit .eqv. fld_radforce_split) then
      active = .true.
      !> Radiation force = kappa*rho/c *Flux = lambda gradE
      call fld_get_fluxlimiter(wct,x,ixi^l,ixo^l,lambda,fld_r,nghostcells)
      do idir = 1,ndim
        call gradient(wct(ixi^s,iw_r_e),ixi^l,ixo^l,idir,grad_e,nghostcells)
        radiation_forcect(ixo^s,idir) = -lambda(ixo^s)*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 
        w(ixo^s,iw_e) = w(ixo^s,iw_e) &
            + qdt*wctprim(ixo^s,iw_mom(idir))*radiation_forcect(ixo^s,idir)
      enddo
    end if
  end subroutine add_fld_rad_force
 
  !> get dt limit for radiation force: NOTE: only uniform cartesian here!
  subroutine fld_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
    double precision :: radiation_force(ixo^s,1:ndim)
    double precision :: dxinv(1:ndim), max_grad
    double precision :: lambda(ixo^s),fld_r(ixo^s),grad_e(ixi^s)
    integer          :: idir
 
    ^d&dxinv(^d)=one/dx^d;
    call fld_get_fluxlimiter(w,x,ixi^l,ixo^l,lambda,fld_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)*grad_e(ixo^s)
      max_grad = maxval(dabs(radiation_force(ixo^s,idir)))
      max_grad = max(max_grad, epsilon(1.0d0))
      dtnew = min(dtnew, courantpar / dsqrt(max_grad * dxinv(idir)))
    end do
  end subroutine fld_radforce_get_dt
 
  !> Sets the opacity in the w-array
  !> by calling mod_opal_opacity
  subroutine fld_get_opacity(w, x, ixI^L, ixO^L, fld_kappa)
    use mod_global_parameters
    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)
    double precision :: temp(ixi^s)
    double precision :: rho0,temp0,n
    integer :: i,j,ix^d, idir
    select case (fld_opacity_law)
      case('const')
        fld_kappa(ixo^s) = fld_kappa0/unit_opacity
      case('thomson')
        fld_kappa(ixo^s) = fld_kappa0/unit_opacity
      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) = n/unit_opacity
        {enddo\ }
      case('special')
        if (.not. associated(usr_special_opacity)) then
          call mpistop("special opacity not defined")
        endif
        call usr_special_opacity(ixi^l, ixo^l, w, x, fld_kappa)
      case default
        call mpistop("Doesn't know opacity law")
      end select
  end subroutine fld_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 fld_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),fld_lambda(ixo^s)
    integer :: idir, ix^d
    integer :: filter
    double precision :: kappa(ixo^s)
    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')
      !> optically thick limit
      fld_lambda(ixo^s) = 1.d0/3.d0
      fld_r(ixo^s) = zero
    case('FreeStream')
      !> optically thin limit
      !> 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,nth)
        normgrad2(ixo^s) = normgrad2(ixo^s)+grad_r_e(ixo^s)**2
      end do
      call fld_get_opacity(w,x,ixi^l,ixo^l,kappa)
      ! Note: possible zeros for uniform setup, to adjust
      if(maxval(normgrad2(ixo^s))==0.0d0)call mpistop("choose other fluxlimiter, FreeStream has issues")
      fld_r(ixo^s) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
      !> Calculate the flux limiter, lambda
      fld_lambda(ixo^s) = one/fld_r(ixo^s)
    case('Pomraning')
      !> 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,nth)
        normgrad2(ixo^s) = normgrad2(ixo^s) + grad_r_e(ixo^s)**2
      end do
      call fld_get_opacity(w,x,ixi^l,ixo^l,kappa)
      fld_r(ixo^s) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s)*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) = (2.d0+fld_r(ixo^s))/(6.d0+3*fld_r(ixo^s)+fld_r(ixo^s)**2.d0)
    case('Minerbo')
      !> 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,nth)
        normgrad2(ixo^s) = normgrad2(ixo^s) + grad_r_e(ixo^s)**2
      end do
      call fld_get_opacity(w, x, ixi^l, ixo^l, kappa)
      fld_r(ixo^s) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s)*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) .lt. 3.d0/2.d0) then
          fld_lambda(ix^d) = 2.d0/(3.d0+dsqrt(9.d0+12.d0*fld_r(ix^d)**2.d0))
        else
          fld_lambda(ix^d) = 1.d0/(1.d0+fld_r(ix^d)+dsqrt(1.d0+2.d0*fld_r(ix^d)))
        endif
      {enddo\}
    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
 
  end subroutine fld_get_fluxlimiter
 
  !> Calculate Radiation Flux
  !> stores radiation flux in w-array
  subroutine fld_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)
    double precision :: cc
    double precision :: grad_r_e(ixi^s)
    double precision :: rad_e(ixi^s)
    double precision :: kappa(ixo^s), lambda(ixo^s), fld_r(ixo^s)
    integer :: ix^d, idir
 
    cc = const_c/unit_velocity
    rad_e(ixi^s) = w(ixi^s, iw_r_e)
    call fld_get_opacity(w, x, ixi^l, ixo^l, kappa)
    call fld_get_fluxlimiter(w, x, ixi^l, ixo^l, lambda, fld_r, nghostcells)
    !> Calculate the Flux using the fld closure relation
    !> F = -c*lambda/(kappa*rho) *grad E
    do idir = 1,ndim
      call gradient(rad_e,ixi^l,ixo^l,idir,grad_r_e,nth)
      rad_flux(ixo^s, idir) = -cc*lambda(ixo^s)/(kappa(ixo^s)*w(ixo^s,iw_rho))*grad_r_e(ixo^s)
    end do
  end subroutine fld_get_radflux
 
  !> Calculate Eddington-tensor
  subroutine fld_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(ixI^S,1:ndim,1:ndim)
    double precision :: tnsr2(ixO^S,1:ndim,1:ndim)
    double precision :: normgrad2(ixO^S), f(ixO^S)
    double precision :: grad_r_e(ixI^S,1:ndim), rad_e(ixI^S)
    double precision :: lambda(ixO^S), fld_R(ixO^S)
    integer :: i,j, idir,jdir
 
    !> 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
    end do
    call fld_get_fluxlimiter(w,x,ixi^l,ixo^l,lambda,fld_r,nth)
    !> Calculate radiation pressure
    !> P = (lambda + lambda^2 R^2)*E
    f(ixo^s) = lambda(ixo^s)+lambda(ixo^s)**two*fld_r(ixo^s)**two
    f(ixo^s) = half*(one-f(ixo^s))+half*(3.d0*f(ixo^s)-one)
    {^ifoned
    eddington_tensor(ixo^s,1,1) = f(ixo^s)
    }
    {^nooned
    do idir = 1,ndim
      eddington_tensor(ixo^s,idir,idir) = half*(one-f(ixo^s))
    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)-one) &
            * 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 fld_get_eddington
 
  !> Calculate Radiation Pressure
  !> Returns Radiation Pressure as tensor
  !> NOTE: w is primitive on entry
  subroutine fld_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(ixi^s,1:ndim,1:ndim)
 
    call fld_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 fld_get_radpress
 
  subroutine fld_implicit_update(dtfactor,qdt,qtC,psa,psb)
    use mod_global_parameters
    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 :: ixO^L,iigrid,igrid
 
    call diffuse_e_rad_mg(dtfactor,qdt,qtc,psa,psb)
 
    ixo^l=ixg^ll^lsubnghostcells;
    !$OMP PARALLEL DO PRIVATE(igrid)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       call energy_interaction(psa(igrid)%w, psa(igrid)%x, ixg^ll, ixo^l, dtfactor, qdt); 
    end do
    !$OMP END PARALLEL DO
  end subroutine fld_implicit_update
 
  !> Energy interaction and photon tiring
  subroutine energy_interaction(w, x, ixI^L, ixO^L, dtfactor, qdt)
    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)    :: dtfactor, qdt
    double precision, intent(inout) :: w(ixI^S,1:nw)
 
    double precision :: wprim(ixI^S,1:nw)
    double precision, dimension(ixO^S) :: a1,a2,a3,c0,c1,kappa
    double precision, dimension(ixI^S) :: e_gas,E_rad,vel,grad_v,nabla_vP
    double precision, dimension(ixI^S,1:ndim,1:ndim) :: div_v,edd
    double precision :: sigma_b,cc
    integer :: i,j,idir,jdir,ix^D
 
    !> Photon tiring
    !> calculate tensor div_v
    !> !$OMP PARALLEL DO
    do idir = 1,ndim
      do jdir = 1,ndim
        vel(ixi^s) = w(ixi^s,iw_mom(jdir))/w(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 fld_get_eddington(w,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)
    }
 
    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) = w(ixo^s,iw_e)-half*sum(w(ixo^s,iw_mom(:))**2,dim=ndim+1)/w(ixo^s,iw_rho)
    if(allocated(iw_mag)) then
      e_gas(ixo^s) = e_gas(ixo^s)-half*sum(w(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/(fld_gamma-1.d0))
    {enddo\}
    e_rad(ixo^s) = w(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 fld_get_opacity(w,x,ixi^l,ixo^l,kappa)
    endif
    !> Coefficients for the polynomial in Moens et al. 2022, eq 37.
    a1(ixo^s) = 4.d0*kappa(ixo^s)*sigma_b*(fld_gamma-one)**4/w(ixo^s,iw_rho)**3*dtfactor*qdt
    a2(ixo^s) = cc*kappa(ixo^s)*w(ixo^s,iw_rho)*dtfactor*qdt
    a3(ixo^s) = nabla_vp(ixo^s)*dtfactor*qdt
    c0(ixo^s) = ((one+a2(ixo^s)+a3(ixo^s))*e_gas(ixo^s)+a2(ixo^s)*e_rad(ixo^s))/a1(ixo^s)/(one+a3(ixo^s))
    c1(ixo^s) = (one+a2(ixo^s)+a3(ixo^s))/a1(ixo^s)/(one+a3(ixo^s))
 
    !> 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 default
        call mpistop('root-method not known')
      end select
    {enddo\}
    e_rad(ixo^s) = (a1(ixo^s)*e_gas(ixo^s)**4.d0+e_rad(ixo^s))/(one+a2(ixo^s)+a3(ixo^s))
    !> 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(w(ixo^s,iw_mom(:))**2,dim=ndim+1)/w(ixo^s,iw_rho)
    if(allocated(iw_mag)) then
      w(ixo^s,iw_e) = w(ixo^s,iw_e)+half*sum(w(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
 
  !> 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_multigrid_coupling
    type(state), target :: psa(max_blocks) !< Advance psa=psb+dtfactor*qdt*F_im(psa)
    type(state), target :: psb(max_blocks)
    double precision, intent(in) :: qdt
    double precision, intent(in) :: qtC
    double precision, intent(in) :: dtfactor
    integer                      :: n
    double precision             :: res, max_residual, lambda
    integer, parameter           :: max_its = 100
    integer                      :: iw_to,iw_from
    integer                      :: iigrid, igrid, id
    integer                      :: nc, lvl, i
    type(tree_node), pointer     :: pnode
    double precision             :: fac, facD
 
    !> Set diffusion timestep, add previous timestep if mg did not converge:
    if(it == 0) dt_diff = 0
    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 = fld_diff_tol
    mg%operator_type = mg_vhelmholtz
    mg%smoother_type = mg_smoother_gs
    call mg_set_methods(mg)
    if(.not. mg%is_allocated) call mpistop("multigrid tree not allocated yet")
    lambda = 1.d0/(dtfactor * dt_diff)
    call vhelmholtz_set_lambda(lambda)
    call update_diffcoeff(psb)
    fac = 1.d0
    facd = 1.d0
    call mg_copy_to_tree(i_diff_mg, mg_iveps, factor=facd, state_from=psb)
    call mg_copy_to_tree(iw_r_e, mg_iphi, factor=fac, state_from=psb)
    call mg_copy_to_tree(iw_r_e, mg_irhs, factor=-1/(dtfactor*dt_diff), state_from=psb)
    if(time_advance)then
      call mg_restrict(mg, mg_iveps)
      call mg_fill_ghost_cells(mg, mg_iveps)
    endif
    call mg_fas_fmg(mg, .true., max_res=res)
    do n = 1, max_its
      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
    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=ixm^ll^ladd1;
    !$OMP PARALLEL DO PRIVATE(igrid)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
       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)*w(hxo^s,i_diff_mg)/(w(ixo^s,i_diff_mg) + w(hxo^s,i_diff_mg))*(w(hxo^s,iw_r_e) - w(ixo^s,iw_r_e))
      divf_j(ixo^s) = w(ixo^s,i_diff_mg)*w(jxo^s,i_diff_mg)/(w(ixo^s,i_diff_mg) + w(jxo^s,i_diff_mg))*(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 fld_get_diffcoef_central(w, wCT, wCTprim, x, ixI^L, ixO^L, primitives_filled)
    use mod_global_parameters
    use mod_geometry
    use mod_usr_methods
    integer, intent(in)          :: ixi^l, ixo^l
    logical, intent(in)          :: primitives_filled
    double precision, intent(inout) :: w(ixi^s,1:nw)
    double precision, intent(in) :: wct(ixi^s,1:nw)
    double precision, intent(in) :: wctprim(ixi^s,1:nw)
    double precision, intent(in) :: x(ixi^s,1:ndim)
    integer :: idir,i,j,ix^d
    double precision :: cc
    double precision :: kappa(ixo^s),lambda(ixo^s),fld_r(ixo^s)
    double precision :: max_d(ixi^s),grad_r_e(ixi^s),rad_e(ixi^s)
 
    cc = const_c/unit_velocity
    if(primitives_filled)then
       ! TODO: pass primitives here
       call fld_get_opacity(wct, x, ixi^l, ixo^l, kappa)
       call fld_get_fluxlimiter(wct, x, ixi^l, ixo^l, lambda, fld_r, nghostcells-1)
       !> calculate diffusion coefficient
       w(ixi^s,i_diff_mg) = zero !> so that w(i_diff_mg) in ghostcells won't accumulate to extreme values
       w(ixo^s,i_diff_mg) = cc*lambda(ixo^s)/(kappa(ixo^s)*wct(ixo^s,iw_rho))
       where(w(ixo^s,i_diff_mg) .lt. 0.d0) w(ixo^s,i_diff_mg) = smalldouble
     else
       ! WARNING: in this case, w, wCT, wCTprim are exactly the same storage location
       call fld_get_opacity(wct, x, ixi^l, ixo^l, kappa)
       call fld_get_fluxlimiter(wct, x, ixi^l, ixo^l, lambda, fld_r, nghostcells-1)
       !> calculate diffusion coefficient
       w(ixi^s,i_diff_mg) = zero !> so that w(i_diff_mg) in ghostcells won't accumulate to extreme values
       w(ixo^s,i_diff_mg) = cc*lambda(ixo^s)/(kappa(ixo^s)*wct(ixo^s,iw_rho))
       where(w(ixo^s,i_diff_mg) .lt. 0.d0) w(ixo^s,i_diff_mg) = smalldouble
     endif
     if(associated(usr_special_diffcoef)) call usr_special_diffcoef(w, wct, x, ixi^l, ixo^l)
 
  end subroutine fld_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=ixm^ll^ladd1;
    !$OMP PARALLEL DO PRIVATE(igrid)
    do iigrid=1,igridstail; igrid=igrids(iigrid);
       ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
        call fld_get_diffcoef_central(psa(igrid)%w, psa(igrid)%w, psa(igrid)%w, psa(igrid)%x, ixg^ll, ixo^l, .false.)
    end do
    !$OMP END PARALLEL DO
  end subroutine update_diffcoeff
 
  !> 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 = 100
    bisect_a = zero
    bisect_b = min(abs(c0/c1),abs(c0)**(1.d0/4.d0))+smalldouble
    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
        print*, 'skip to bisection algorithm'
        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_fld