advance_species_transport Subroutine

public subroutine advance_species_transport(mesh, flow, bc, params, fields, species, transport)

Performs one explicit Euler step for non-reacting species transport.

This routine evaluates: - Constant-density advection: fields%face_flux [m^3/s] advects directly. - Variable-density advection: fields%mass_flux [kg/s] advances conservative , then divides by the active density. - Diffusive fluxes: constant-density mode uses ; variable-density mode uses . - Correction velocity: subtracts at each face so the corrected diffusive species flux has zero net mass flux. - Bounding & Normalization: Clamps and enforces locally. Reactions are not applied here; enable_reactions is rejected by the supported variable-density path.

Arguments

Type IntentOptional Attributes Name
type(mesh_t), intent(in) :: mesh

The computational mesh.
type(flow_mpi_t), intent(inout) :: flow

MPI decomposition data for synchronization.
type(bc_set_t), intent(in) :: bc

Boundary condition settings.
type(case_params_t), intent(in) :: params

Simulation parameters (dt, etc).
type(flow_fields_t), intent(in) :: fields

Flow field (velocity/face fluxes).
type(species_fields_t), intent(inout) :: species

Mass fraction fields to update.
type(transport_properties_t), intent(in) :: transport

Physical properties (diffusivities).

Source Code

   subroutine advance_species_transport(mesh, flow, bc, params, fields, species, transport)
      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(species_fields_t), intent(inout) :: species
      type(transport_properties_t), intent(in) :: transport

      real(rk), allocatable :: dY(:), diff_flux(:), adv_flux(:), Y_face_lin(:), rhs_current(:)
      real(rk), allocatable :: gradY(:,:,:)
      real(rk) :: flux, face_area, dist, sum_diff_flux
      real(rk) :: rho_face, rho_cell, rho_old_cell
      real(rk) :: Y_cell, Y_other, Y_face
      real(rk) :: diff_face, rhs_update
      real(rk) :: c_curr, c_old
      integer :: c, f, fid, neighbor
      integer :: k
      logical :: is_dirichlet, use_variable_density_species, use_ab2

      if (species%nspecies <= 0) return

      use_ab2 = trim(lowercase(params%species_time_scheme)) == 'ab2' .or. &
                trim(lowercase(params%species_time_scheme)) == 'adams_bashforth2' .or. &
                trim(lowercase(params%species_time_scheme)) == 'adams-bashforth2'
      c_curr = 1.5_rk
      c_old = -0.5_rk
      if (use_ab2 .and. species%rhs_old_valid .and. params%dt_old > tiny_safe) then
         c_curr = 1.0_rk + 0.5_rk * (params%dt / params%dt_old)
         c_old = -0.5_rk * (params%dt / params%dt_old)
      end if
      use_variable_density_species = params%enable_variable_density
      if (use_variable_density_species) then
         if (.not. allocated(fields%mass_flux)) then
            call fatal_error('species', 'variable-density species transport requires fields%mass_flux')
         end if
         if (.not. allocated(transport%rho)) then
            call fatal_error('species', 'variable-density species transport requires transport%rho')
         end if
         if (.not. allocated(transport%rho_old)) then
            call fatal_error('species', 'variable-density species transport requires transport%rho_old')
         end if
      end if

      allocate(dY(species%nspecies))
      allocate(diff_flux(species%nspecies))
      allocate(adv_flux(species%nspecies))
      allocate(Y_face_lin(species%nspecies))
      allocate(rhs_current(species%nspecies))
      allocate(gradY(3, mesh%ncells, species%nspecies))

      species%Y_old = species%Y
      do k = 1, species%nspecies
         call compute_scalar_gradients(mesh, species%Y_old(k, :), gradY(:, :, k))
      end do

      ! Iterate through MPI-owned cells
      do c = flow%first_cell, flow%last_cell
         dY = zero

         do f = 1, mesh%ncell_faces(c)
            fid = mesh%cell_faces(f,c)
            if (mesh%faces(fid)%owner == c) then
               if (use_variable_density_species) then
                  flux = fields%mass_flux(fid)
               else
                  flux = fields%face_flux(fid)
               end if
            else
               if (use_variable_density_species) then
                  flux = -fields%mass_flux(fid)
               else
                  flux = -fields%face_flux(fid)
               end if
            end if
            neighbor = face_effective_neighbor(mesh, bc, fid, c)

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

            sum_diff_flux = zero
            do k = 1, species%nspecies
               Y_cell = species%Y_old(k, c)
               
               if (neighbor == 0) then
                  call boundary_species(mesh, bc, fid, k, Y_cell, Y_other, is_dirichlet)
               else
                  Y_other = species%Y_old(k, neighbor)
                  is_dirichlet = .true.
               end if

               ! 1. Advective flux.  The selected scalar scheme can be
               ! first-order upwind, central, bounded central, or compact
               ! limited-linear/MUSCL reconstruction.  flux is oriented
               ! outward from the current cell.
               Y_face = scalar_face_value(mesh, species%Y_old(k, :), gradY(:, :, k), &
                                          params%species_convection_scheme, params%scalar_limiter, &
                                          c, neighbor, fid, flux, Y_other, is_dirichlet)
               adv_flux(k) = -flux * Y_face

               ! 2. Diffusive flux using central difference
               diff_flux(k) = zero
               if (neighbor /= 0 .or. is_dirichlet) then
                  if (neighbor == 0) then
                     diff_face = transport%diffusivity(k, c)
                  else
                     diff_face = 0.5_rk * (transport%diffusivity(k, c) + transport%diffusivity(k, neighbor))
                  end if
                  if (use_variable_density_species) then
                     if (neighbor == 0) then
                        rho_face = transport%rho(c)
                     else
                        rho_face = 0.5_rk * (transport%rho(c) + transport%rho(neighbor))
                     end if
                     diff_flux(k) = rho_face * diff_face * (Y_other - Y_cell) / dist * face_area
                  else
                     diff_flux(k) = diff_face * (Y_other - Y_cell) / dist * face_area
                  end if
               end if
               
               sum_diff_flux = sum_diff_flux + diff_flux(k)
               Y_face_lin(k) = 0.5_rk * (Y_cell + Y_other)
            end do

            ! 3. Apply Correction Velocity to ensure mass conservation: sum(j_k) = 0
            do k = 1, species%nspecies
               dY(k) = dY(k) + adv_flux(k) + (diff_flux(k) - Y_face_lin(k) * sum_diff_flux)
            end do
         end do

         ! Explicit timestep update. Constant-density mode updates Y directly.
         ! Variable-density mode updates rho*Y using mass fluxes, then divides
         ! by the active cell density for this split substep.
         rho_cell = one
         rho_old_cell = one
         if (use_variable_density_species) then
            rho_cell = max(transport%rho(c), tiny_safe)
            rho_old_cell = max(transport%rho_old(c), tiny_safe)
         end if

         rhs_current = dY / mesh%cells(c)%volume

         do k = 1, species%nspecies
            if (use_ab2 .and. species%rhs_old_valid) then
               rhs_update = c_curr * rhs_current(k) + c_old * species%rhs_old(k, c)
            else
               rhs_update = rhs_current(k)
            end if

            if (use_variable_density_species) then
               species%Y(k,c) = (rho_old_cell * species%Y_old(k,c) + params%dt * rhs_update) / rho_cell
            else
               species%Y(k,c) = species%Y_old(k,c) + params%dt * rhs_update
            end if
            
         end do

         species%rhs_old(:, c) = rhs_current

         ! Keep the local composition on the bounded simplex.  This still is a
         ! fallback repair for explicit transport overshoots, but it is less
         ! destructive than scaling every species after clipping: the sum defect
         ! is preferentially put into/removes from the dominant bath component,
         ! preserving trace-species ratios whenever possible.
         call repair_cell_mass_fractions(species%Y(:, c), species%nspecies)
      end do

      deallocate(dY)
      deallocate(diff_flux)
      deallocate(adv_flux)
      deallocate(Y_face_lin)
      deallocate(rhs_current)
      deallocate(gradY)

      species%rhs_old_valid = .true.

      ! Synchronize updated species ghosts for the next transport/property step.
      call flow_exchange_cell_matrix(flow, species%Y)

   end subroutine advance_species_transport