subroutine radiation_update_source(mesh, bc, flow, rad, context, params, energy, species_Y, step, time, dt)
type(mesh_t), intent(in) :: mesh
type(bc_set_t), intent(in) :: bc
type(flow_mpi_t), intent(inout) :: flow
type(radiation_mpi_t), intent(in) :: rad
type(radiation_context_t), intent(inout) :: context
type(case_params_t), intent(in) :: params
type(energy_fields_t), intent(inout) :: energy
real(rk), intent(in), optional :: species_Y(:,:)
integer, intent(in) :: step
real(rk), intent(in) :: time, dt
type(radiation_state_t) :: state
type(radiation_source_t) :: partial, total
real(rk) :: t_start, t_after_gather, t_after_compute, t_after_reduce
real(rk) :: gather_time, compute_time, reduce_time
real(rk) :: debug_T_diff, debug_Y_diff
real(rk) :: absorption
real(rk), parameter :: sigma = 5.670374419e-8_rk
integer :: c, ierr
if (.not. context%enabled) return
if (.not. params%enable_radiation) return
if (.not. allocated(energy%qrad)) return
if (params%radiation_update_interval <= 0) return
if (mod(step-1, params%radiation_update_interval) /= 0 .and. step /= 1) return
call profiler_start('Radiation_Source_Update')
t_start = real(MPI_Wtime(), rk)
call allocate_radiation_state(state, mesh%ncells, context%n_species, context%n_scalars)
call allocate_radiation_source(partial, mesh%ncells)
call allocate_radiation_source(total, mesh%ncells)
call profiler_start('Radiation_State_Gather')
call prepare_radiation_state(mesh, flow, context, params, energy, species_Y, step, time, dt, state, &
debug_T_diff, debug_Y_diff)
call profiler_stop('Radiation_State_Gather')
t_after_gather = real(MPI_Wtime(), rk)
call profiler_start('Radiation_Model_Compute')
select case (trim(context%model))
case ('none')
call radiation_compute_none(mesh, context, state, partial)
case ('spectral_test')
call radiation_compute_spectral_test(mesh, context, state, partial, params%radiation_source_scale)
case ('external')
call radiation_compute_external(mesh, context, state, partial)
case ('p1')
call radiation_compute_p1(mesh, bc, context, state, partial, params)
case ('dom')
call radiation_compute_dom(mesh, bc, context, state, partial, params)
case default
call fatal_error('radiation', 'unknown radiation_source_model: '//trim(context%model))
end select
call profiler_stop('Radiation_Model_Compute')
t_after_compute = real(MPI_Wtime(), rk)
! Retrieve the unreduced partial absorption term (kappa*G or kappa*G_local) from
! each wavenumber spectral segment.
total%qrad = zero
call profiler_start('Radiation_Source_Reduce')
call MPI_Allreduce(partial%qrad, total%qrad, mesh%ncells, MPI_DOUBLE_PRECISION, MPI_SUM, rad%comm, ierr)
if (ierr /= MPI_SUCCESS) call fatal_error('radiation', 'MPI failure reducing radiation qrad')
call profiler_stop('Radiation_Source_Reduce')
! [X-1 Contract]: In radiation model modules, source%qrad stores only the unreduced
! partial absorption contribution kappa*G (or kappa*G_local).
! The radiation driver is responsible for performing the MPI reduction above and then
! subtracting the isotropic emission contribution 4*kappa*sigma*T^4 here locally on
! every rank to form the final net radiative gas energy source:
! q_rad = kappa*G - 4*kappa*sigma*T^4.
if (trim(context%model) == 'dom' .or. trim(context%model) == 'p1') then
absorption = max(params%radiation_absorption_coeff, 1.0e-6_rk)
do c = 1, mesh%ncells
total%qrad(c) = total%qrad(c) - absorption * 4.0_rk * sigma * (state%temperature(c)**4)
end do
end if
energy%qrad = total%qrad
t_after_reduce = real(MPI_Wtime(), rk)
gather_time = t_after_gather - t_start
compute_time = t_after_compute - t_after_gather
reduce_time = t_after_reduce - t_after_compute
call write_radiation_diagnostics(mesh, rad, context, state, total, step, time, dt, &
gather_time, compute_time, reduce_time, debug_T_diff, debug_Y_diff)
call free_radiation_source(total)
call free_radiation_source(partial)
call free_radiation_state(state)
call profiler_stop('Radiation_Source_Update')
end subroutine radiation_update_source