advance_energy_transport Subroutine

public subroutine advance_energy_transport(mesh, flow, bc, params, fields, energy, species_Y)

Advance transported sensible enthalpy with constant flow density.

The explicit finite-volume update is:

V dh/dt = - sum_f F_f h_f + (1/rho) sum_f lambda A_f (T_nb - T_c)/d_f + (qrad/rho) V

fields%face_flux is volumetric flux. It is re-oriented outward from the currently updated cell before applying upwind h advection. Cantera thermo may provide cp and lambda, but the projection/momentum density remains params%rho.

Arguments

Type IntentOptional Attributes Name
type(mesh_t), intent(in) :: mesh
type(flow_mpi_t), intent(inout) :: flow
type(bc_set_t), intent(in) :: bc
type(case_params_t), intent(in) :: params
type(flow_fields_t), intent(in) :: fields
type(energy_fields_t), intent(inout) :: energy
real(kind=rk), intent(in), optional :: species_Y(:,:)

Calls

proc~~advance_energy_transport~~CallsGraph proc~advance_energy_transport mod_energy::advance_energy_transport proc~boundary_enthalpy_from_temperature mod_energy::boundary_enthalpy_from_temperature proc~advance_energy_transport->proc~boundary_enthalpy_from_temperature proc~boundary_temperature mod_bc::boundary_temperature proc~advance_energy_transport->proc~boundary_temperature proc~build_boundary_thermo_y mod_energy::build_boundary_thermo_Y proc~advance_energy_transport->proc~build_boundary_thermo_y proc~energy_face_normal_distance mod_energy::energy_face_normal_distance proc~advance_energy_transport->proc~energy_face_normal_distance proc~face_effective_neighbor mod_bc::face_effective_neighbor proc~advance_energy_transport->proc~face_effective_neighbor proc~fatal_error mod_kinds::fatal_error proc~advance_energy_transport->proc~fatal_error proc~flow_exchange_cell_scalar mod_mpi_flow::flow_exchange_cell_scalar proc~advance_energy_transport->proc~flow_exchange_cell_scalar proc~profiler_start mod_profiler::profiler_start proc~advance_energy_transport->proc~profiler_start proc~profiler_stop mod_profiler::profiler_stop proc~advance_energy_transport->proc~profiler_stop proc~recover_temperature_and_update_thermo_cantera mod_energy::recover_temperature_and_update_thermo_cantera proc~advance_energy_transport->proc~recover_temperature_and_update_thermo_cantera proc~recover_temperature_constant_cp mod_energy::recover_temperature_constant_cp proc~advance_energy_transport->proc~recover_temperature_constant_cp proc~enthalpy_from_temperature_cantera_point mod_energy::enthalpy_from_temperature_cantera_point proc~boundary_enthalpy_from_temperature->proc~enthalpy_from_temperature_cantera_point proc~enthalpy_from_temperature_value mod_energy::enthalpy_from_temperature_value proc~boundary_enthalpy_from_temperature->proc~enthalpy_from_temperature_value proc~build_boundary_thermo_y->proc~fatal_error proc~boundary_species mod_bc::boundary_species proc~build_boundary_thermo_y->proc~boundary_species proc~energy_face_normal_distance->proc~fatal_error proc~energy_outward_normal mod_energy::energy_outward_normal proc~energy_face_normal_distance->proc~energy_outward_normal proc~patch_type_for_face mod_bc::patch_type_for_face proc~energy_face_normal_distance->proc~patch_type_for_face proc~is_periodic_face mod_bc::is_periodic_face proc~face_effective_neighbor->proc~is_periodic_face proc~flow_exchange_cell_scalar->proc~profiler_start proc~flow_exchange_cell_scalar->proc~profiler_stop mpi_irecv mpi_irecv proc~flow_exchange_cell_scalar->mpi_irecv mpi_isend mpi_isend proc~flow_exchange_cell_scalar->mpi_isend mpi_waitall mpi_waitall proc~flow_exchange_cell_scalar->mpi_waitall proc~check_mpi~2 mod_mpi_flow::check_mpi proc~flow_exchange_cell_scalar->proc~check_mpi~2 mpi_wtime mpi_wtime proc~profiler_start->mpi_wtime proc~find_or_create_timer mod_profiler::find_or_create_timer proc~profiler_start->proc~find_or_create_timer proc~profiler_stop->mpi_wtime proc~profiler_stop->proc~find_or_create_timer proc~record_edge mod_profiler::record_edge proc~profiler_stop->proc~record_edge proc~recover_temperature_and_update_thermo_cantera->proc~fatal_error interface~cantera_recover_temperature_and_update_thermo_c mod_energy::cantera_recover_temperature_and_update_thermo_c proc~recover_temperature_and_update_thermo_cantera->interface~cantera_recover_temperature_and_update_thermo_c proc~build_c_species_names mod_energy::build_c_species_names proc~recover_temperature_and_update_thermo_cantera->proc~build_c_species_names proc~build_thermo_y mod_energy::build_thermo_Y proc~recover_temperature_and_update_thermo_cantera->proc~build_thermo_y proc~recover_temperature_constant_cp->proc~fatal_error proc~build_thermo_y->proc~fatal_error proc~check_mpi~2->proc~fatal_error proc~enthalpy_from_temperature_cantera_point->proc~fatal_error proc~enthalpy_from_temperature_cantera_point->proc~build_c_species_names interface~cantera_update_thermo_c mod_energy::cantera_update_thermo_c proc~enthalpy_from_temperature_cantera_point->interface~cantera_update_thermo_c proc~is_periodic_face->proc~patch_type_for_face

Called by

proc~~advance_energy_transport~~CalledByGraph proc~advance_energy_transport mod_energy::advance_energy_transport program~lowmach_react_hex lowmach_react_hex program~lowmach_react_hex->proc~advance_energy_transport

Source Code

   subroutine advance_energy_transport(mesh, flow, bc, params, fields, energy, species_Y)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(inout) :: flow
      type(bc_set_t), intent(in) :: bc
      type(case_params_t), intent(in) :: params
      type(flow_fields_t), intent(in) :: fields
      type(energy_fields_t), intent(inout) :: energy
      real(rk), intent(in), optional :: species_Y(:,:)

      integer :: c, lf, fid, neighbor
      real(rk) :: flux_out
      real(rk) :: face_area, dist
      real(rk) :: h_cell, h_other, h_face
      real(rk) :: T_cell, T_other
      real(rk) :: rhs_h, diff_term, lambda_face
      logical :: is_dirichlet, do_diffusion
      real(rk), allocatable :: Y_point(:)

      if (.not. params%enable_energy) return

      if (.not. energy%initialized) then
         call fatal_error('energy', 'advance requested before energy initialization')
      end if

      if (params%rho <= tiny_safe) call fatal_error('energy', 'rho must be positive')
      if (params%energy_cp <= tiny_safe) call fatal_error('energy', 'energy_cp must be positive')
      if (params%energy_lambda < zero) call fatal_error('energy', 'energy_lambda must be non-negative')

      ! Option A: preserve transported enthalpy across species updates.
      ! Species transport may have changed Y before this energy step.  Do not
      ! recompute h from the old temperature and the new composition.  Instead,
      ! keep h as the transported state, recover T(h,Y,p0), then refresh the
      ! thermo/transport properties used by the enthalpy equation.
      call profiler_start('Energy_Exchange_H')
      call flow_exchange_cell_scalar(flow, energy%h)
      call profiler_stop('Energy_Exchange_H')

      energy%h_old = energy%h

      if (params%enable_cantera_thermo) then
         ! A pre-flux thermo sync is required only when species transport may
         ! have changed the composition since the previous post-flux sync.
         ! Without species, current T/cp/lambda/rho_thermo are already valid
         ! from initialization or the previous energy step's post-sync.
         if (params%enable_species .and. present(species_Y) .and. params%nspecies > 0) then
            call profiler_start('Energy_Cantera_PreSync')
            call recover_temperature_and_update_thermo_cantera(mesh, params, energy, species_Y)
            call profiler_stop('Energy_Cantera_PreSync')

            energy%h = energy%h_old
         end if
      end if

      ! Ensure off-rank temperature/property values are current before fluxes.
      call profiler_start('Energy_PreFlux_Exchange')
      call flow_exchange_cell_scalar(flow, energy%T)
      if (allocated(energy%lambda)) call flow_exchange_cell_scalar(flow, energy%lambda)
      call profiler_stop('Energy_PreFlux_Exchange')

      if (params%enable_cantera_thermo .and. params%enable_species .and. &
          present(species_Y) .and. params%nspecies > 0) then
         allocate(Y_point(params%nspecies))
      end if
      if (allocated(energy%T_old)) energy%T_old = energy%T

      call profiler_start('Energy_Flux_Update')
      do c = flow%first_cell, flow%last_cell
         rhs_h = zero
         h_cell = energy%h_old(c)
         T_cell = energy%T(c)

         do lf = 1, mesh%ncell_faces(c)
            fid = mesh%cell_faces(lf, c)

            ! fields%face_flux is oriented owner -> neighbor. Reorient outward
            ! from the current cell.
            if (mesh%faces(fid)%owner == c) then
               flux_out = fields%face_flux(fid)
            else
               flux_out = -fields%face_flux(fid)
            end if

            neighbor = face_effective_neighbor(mesh, bc, fid, c)
            face_area = mesh%faces(fid)%area

            if (neighbor > 0) then
               h_other = energy%h_old(neighbor)
               T_other = energy%T(neighbor)
               do_diffusion = .true.
            else
               call boundary_temperature(mesh, bc, fid, T_cell, T_other, is_dirichlet)
               if (is_dirichlet) then
                  if (allocated(Y_point)) then
                     call build_boundary_thermo_Y(mesh, bc, params, fid, c, species_Y, Y_point)
                     h_other = boundary_enthalpy_from_temperature(mesh, params, T_other, Y_point)
                  else
                     h_other = boundary_enthalpy_from_temperature(mesh, params, T_other)
                  end if
                  do_diffusion = .true.
               else
                  h_other = h_cell
                  T_other = T_cell
                  do_diffusion = .false.
               end if
            end if

            ! Upwind advection of h using outward volumetric flux.
            if (flux_out >= zero) then
               h_face = h_cell
            else
               h_face = h_other
            end if
            rhs_h = rhs_h - flux_out * h_face

            ! Fourier heat conduction uses grad(T), not grad(h).
            if (do_diffusion .and. (params%enable_cantera_thermo .or. params%energy_lambda > zero)) then
               dist = energy_face_normal_distance(mesh, bc, fid, c, neighbor)

               if (allocated(energy%lambda)) then
                  if (neighbor > 0) then
                     lambda_face = 0.5_rk * (energy%lambda(c) + energy%lambda(neighbor))
                  else
                     lambda_face = energy%lambda(c)
                  end if
               else
                  lambda_face = params%energy_lambda
               end if

               if (lambda_face > zero) then
                  diff_term = (lambda_face / params%rho) * &
                              (T_other - T_cell) / dist * face_area
                  rhs_h = rhs_h + diff_term
               end if
            end if
         end do

         energy%h(c) = energy%h_old(c) + params%dt * &
                       (rhs_h / mesh%cells(c)%volume + energy%qrad(c) / params%rho)
      end do
      call profiler_stop('Energy_Flux_Update')

      if (params%enable_cantera_thermo) then
         ! Protect transported h from property-refresh roundoff: recover T from
         ! the newly transported h, refresh cp/lambda/rho_thermo, then restore h.
         energy%h_old = energy%h

         call profiler_start('Energy_Cantera_PostSync')
         call recover_temperature_and_update_thermo_cantera(mesh, params, energy, species_Y)
         call profiler_stop('Energy_Cantera_PostSync')

         energy%h = energy%h_old
      else
         call profiler_start('Energy_ConstantCp_RecoverT')
         call recover_temperature_constant_cp(params, energy)
         call profiler_stop('Energy_ConstantCp_RecoverT')
      end if

      ! Synchronize updated owned cells for output and the next step.
      call profiler_start('Energy_Final_Exchange')
      call flow_exchange_cell_scalar(flow, energy%h)
      call flow_exchange_cell_scalar(flow, energy%T)
      if (allocated(energy%lambda)) call flow_exchange_cell_scalar(flow, energy%lambda)
      call profiler_stop('Energy_Final_Exchange')
   end subroutine advance_energy_transport