Advance transported sensible enthalpy one explicit finite-volume step.
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 in .
fields%mass_flux is density-weighted flux in .
Both are stored owner-to-neighbor and re-oriented outward from the
currently updated cell before upwind advection.
Constant-density mode updates h with params%rho. Variable-density
mode uses a conservative rho*h update with fields%mass_flux,
transport%rho_old, and current transport%rho, then stores the
transported specific enthalpy back in energy%h.
| Type | Intent | Optional | 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(:,:) | ||
| type(transport_properties_t), | intent(in), | optional | :: | transport | ||
| integer, | intent(in), | optional | :: | step |
subroutine advance_energy_transport(mesh, flow, bc, params, fields, energy, species_Y, transport, step) 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(:,:) type(transport_properties_t), intent(in), optional :: transport integer, intent(in), optional :: step integer :: c, lf, fid, neighbor, k, ierr real(rk) :: flux_out real(rk) :: face_area, dist, vol real(rk) :: h_cell, h_other, h_face real(rk) :: T_cell, T_other real(rk) :: rhs_h, rhs_h_current, rhs_h_update, diff_term, lambda_face real(rk) :: c_curr, c_old real(rk) :: conductive_boundary_flux_out_local real(rk) :: advective_boundary_rho_h_flux_out_local real(rk) :: qrad_integral_step_local real(rk) :: rho_species_hdiff_integral_step_local real(rk) :: rho_diag real(rk) :: energy_update_delta_rate_local, energy_update_rhs_integral_local real(rk) :: energy_update_delta_rate_global, energy_update_rhs_integral_global real(rk) :: energy_update_balance_global, energy_update_denom real(rk) :: update_send(3), update_recv(3) real(rk) :: operator_consistent_rho_h_integral_local real(rk) :: operator_consistent_rho_h_integral_global real(rk) :: rho_cell, rho_old_cell, rho_face real(rk) :: sum_diff_flux, diff_face, species_hdiff_rhs real(rk) :: Y_cell, Y_other logical :: is_dirichlet, do_diffusion, use_species_hdiff, has_species_diffusion_face logical :: use_variable_density_energy, use_ab2 real(rk), allocatable :: Y_point(:) ! Explicitly null-initialized pointer scratch avoids false gfortran ! release-build descriptor warnings on the conditional h_k diffusion path. real(rk), pointer :: hk_species(:,:) => null(), hk_boundary(:) => null() real(rk), allocatable :: sp_diff_flux(:), sp_Y_face_lin(:) real(rk), allocatable :: grad_h(:,:) 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') use_ab2 = trim(params%energy_time_scheme) == 'ab2' .or. & trim(params%energy_time_scheme) == 'adams_bashforth2' .or. & trim(params%energy_time_scheme) == 'adams-bashforth2' c_curr = 1.5_rk c_old = -0.5_rk if (use_ab2 .and. energy%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_energy = params%enable_variable_density if (use_variable_density_energy) then if (.not. present(transport)) then call fatal_error('energy', 'variable-density energy transport requires transport properties') end if if (.not. allocated(transport%rho)) then call fatal_error('energy', 'variable-density energy transport requires transport%rho') end if if (.not. allocated(transport%rho_old)) then call fatal_error('energy', 'variable-density energy transport requires transport%rho_old') end if if (.not. allocated(fields%mass_flux)) then call fatal_error('energy', 'variable-density energy transport requires fields%mass_flux') end if end if use_species_hdiff = params%enable_species_enthalpy_diffusion if (use_species_hdiff) then if (.not. present(species_Y)) then call fatal_error('energy', 'species enthalpy diffusion requires species_Y') end if if (.not. present(transport)) then call fatal_error('energy', 'species enthalpy diffusion requires transport properties') end if if (.not. allocated(transport%diffusivity)) then call fatal_error('energy', 'species enthalpy diffusion requires transport diffusivity') end if if (.not. params%enable_cantera_thermo) then call fatal_error('energy', 'species enthalpy diffusion requires Cantera thermo') end if if (params%nspecies <= 0) then call fatal_error('energy', 'species enthalpy diffusion requires nspecies > 0') end if if (size(species_Y, 1) < params%nspecies .or. size(species_Y, 2) < mesh%ncells) then call fatal_error('energy', 'species enthalpy diffusion species_Y has incompatible shape') end if if (size(transport%diffusivity, 1) < params%nspecies .or. & size(transport%diffusivity, 2) < mesh%ncells) then call fatal_error('energy', 'species enthalpy diffusion diffusivity has incompatible shape') end if end if ! 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 .and. params%enable_species .and. & present(species_Y) .and. params%nspecies > 0) then allocate(Y_point(params%nspecies)) end if if (use_species_hdiff) then if (.not. allocated(Y_point)) allocate(Y_point(params%nspecies)) allocate(hk_species(params%nspecies, mesh%ncells)) allocate(hk_boundary(params%nspecies)) allocate(sp_diff_flux(params%nspecies)) allocate(sp_Y_face_lin(params%nspecies)) hk_species = zero hk_boundary = zero sp_diff_flux = zero sp_Y_face_lin = zero end if 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 write_variable_density_energy_debug(mesh, flow, params, fields, energy, & 'pre_cantera_presync', species_Y, transport, step) if (use_species_hdiff) then call profiler_start('Energy_Cantera_PreSyncSpeciesH') call recover_temperature_update_thermo_and_species_h_cantera_owned(mesh, flow, params, energy, species_Y, hk_species) call profiler_stop('Energy_Cantera_PreSyncSpeciesH') else call profiler_start('Energy_Cantera_PreSync') call recover_temperature_and_update_thermo_cantera_owned(mesh, flow, params, energy, species_Y) call profiler_stop('Energy_Cantera_PreSync') end if 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') if (allocated(energy%lambda)) then call flow_exchange_cell_scalars(flow, energy%T, energy%lambda) else call flow_exchange_cell_scalar(flow, energy%T) end if call profiler_stop('Energy_PreFlux_Exchange') if (allocated(energy%T_old)) energy%T_old = energy%T allocate(grad_h(3, mesh%ncells)) call compute_scalar_gradients(mesh, energy%h_old, grad_h) if (use_species_hdiff .and. .not. associated(hk_species)) then call fatal_error('energy', 'species enthalpy diffusion scratch was not initialized') end if if (allocated(energy%species_enthalpy_diffusion)) energy%species_enthalpy_diffusion = zero conductive_boundary_flux_out_local = zero energy_update_delta_rate_local = zero energy_update_rhs_integral_local = zero operator_consistent_rho_h_integral_local = zero if (allocated(energy%operator_consistent_rho_h)) energy%operator_consistent_rho_h = zero energy%operator_consistent_rho_h_available = 0 advective_boundary_rho_h_flux_out_local = zero qrad_integral_step_local = zero rho_species_hdiff_integral_step_local = zero 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` and `fields%mass_flux` are oriented ! owner -> neighbor. Reorient outward from the current cell. ! Variable-density energy uses mass flux for conservative rho*h ! advection; constant-density energy uses volumetric flux. if (mesh%faces(fid)%owner == c) then if (use_variable_density_energy) then flux_out = fields%mass_flux(fid) else flux_out = fields%face_flux(fid) end if else if (use_variable_density_energy) then flux_out = -fields%mass_flux(fid) else flux_out = -fields%face_flux(fid) end if end if neighbor = face_effective_neighbor(mesh, bc, fid, c) face_area = mesh%faces(fid)%area is_dirichlet = .false. 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 ! Advective h flux using the selected scalar scheme. The flux ! is outward from the current cell and is volumetric in the ! constant-density path, mass-weighted in variable-density mode. h_face = scalar_face_value(mesh, energy%h_old, grad_h, & params%energy_convection_scheme, params%scalar_limiter, & c, neighbor, fid, flux_out, h_other, is_dirichlet) rhs_h = rhs_h - flux_out * h_face if (neighbor <= 0) then if (use_variable_density_energy) then advective_boundary_rho_h_flux_out_local = advective_boundary_rho_h_flux_out_local + flux_out * h_face else advective_boundary_rho_h_flux_out_local = advective_boundary_rho_h_flux_out_local + params%rho * flux_out * h_face end if end if ! 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 if (use_variable_density_energy) then diff_term = lambda_face * (T_other - T_cell) / dist * face_area else diff_term = (lambda_face / params%rho) * & (T_other - T_cell) / dist * face_area end if rhs_h = rhs_h + diff_term if (neighbor <= 0) then if (use_variable_density_energy) then conductive_boundary_flux_out_local = conductive_boundary_flux_out_local - diff_term else conductive_boundary_flux_out_local = conductive_boundary_flux_out_local - params%rho * diff_term end if end if end if end if ! Species-enthalpy diffusion correction: ! corrected_diff_flux(k) matches mod_species and equals -J_k,out*A/rho. ! Therefore sum_k h_k,face * corrected_diff_flux(k) is added to the ! sensible-enthalpy RHS for this cell. if (use_species_hdiff) then dist = energy_face_normal_distance(mesh, bc, fid, c, neighbor) sum_diff_flux = zero has_species_diffusion_face = .false. do k = 1, params%nspecies Y_cell = species_Y(k, c) if (neighbor > 0) then Y_other = species_Y(k, neighbor) is_dirichlet = .true. else call boundary_species(mesh, bc, fid, k, Y_cell, Y_other, is_dirichlet) end if sp_diff_flux(k) = zero sp_Y_face_lin(k) = 0.5_rk * (Y_cell + Y_other) if (neighbor /= 0 .or. is_dirichlet) then has_species_diffusion_face = .true. 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_energy) then if (neighbor == 0) then rho_face = transport%rho(c) else rho_face = 0.5_rk * (transport%rho(c) + transport%rho(neighbor)) end if sp_diff_flux(k) = rho_face * diff_face * (Y_other - Y_cell) / dist * face_area else sp_diff_flux(k) = diff_face * (Y_other - Y_cell) / dist * face_area end if end if sum_diff_flux = sum_diff_flux + sp_diff_flux(k) end do species_hdiff_rhs = zero if (has_species_diffusion_face) then if (neighbor > 0) then do k = 1, params%nspecies species_hdiff_rhs = species_hdiff_rhs + & 0.5_rk * (hk_species(k, c) + hk_species(k, neighbor)) * & (sp_diff_flux(k) - sp_Y_face_lin(k) * sum_diff_flux) end do else call build_boundary_thermo_Y(mesh, bc, params, fid, c, species_Y, Y_point) call species_sensible_enthalpies_from_temperature_cantera_point(params, T_other, hk_boundary, Y_point) do k = 1, params%nspecies species_hdiff_rhs = species_hdiff_rhs + & 0.5_rk * (hk_species(k, c) + hk_boundary(k)) * & (sp_diff_flux(k) - sp_Y_face_lin(k) * sum_diff_flux) end do end if end if rhs_h = rhs_h + species_hdiff_rhs if (allocated(energy%species_enthalpy_diffusion)) then if (use_variable_density_energy) then rho_cell = max(transport%rho(c), tiny_safe) energy%species_enthalpy_diffusion(c) = energy%species_enthalpy_diffusion(c) + & species_hdiff_rhs / & (rho_cell * mesh%cells(c)%volume) else energy%species_enthalpy_diffusion(c) = energy%species_enthalpy_diffusion(c) + & species_hdiff_rhs / mesh%cells(c)%volume end if end if end if end do if (use_variable_density_energy) then rho_cell = max(transport%rho(c), tiny_safe) rho_old_cell = max(transport%rho_old(c), tiny_safe) rhs_h_current = rhs_h / mesh%cells(c)%volume + energy%qrad(c) if (use_ab2 .and. energy%rhs_old_valid) then rhs_h_update = c_curr * rhs_h_current + c_old * energy%rhs_old(c) else rhs_h_update = rhs_h_current end if energy%h(c) = (rho_old_cell * energy%h_old(c) + params%dt * rhs_h_update) / rho_cell energy_update_delta_rate_local = energy_update_delta_rate_local + & ((rho_cell * energy%h(c) - rho_old_cell * energy%h_old(c)) * mesh%cells(c)%volume) / params%dt energy_update_rhs_integral_local = energy_update_rhs_integral_local + & rhs_h_update * mesh%cells(c)%volume operator_consistent_rho_h_integral_local = operator_consistent_rho_h_integral_local + & rho_cell * energy%h(c) * mesh%cells(c)%volume if (allocated(energy%operator_consistent_rho_h)) then energy%operator_consistent_rho_h(c) = rho_cell * energy%h(c) end if else rhs_h_current = rhs_h / mesh%cells(c)%volume + energy%qrad(c) / params%rho if (use_ab2 .and. energy%rhs_old_valid) then rhs_h_update = c_curr * rhs_h_current + c_old * energy%rhs_old(c) else rhs_h_update = rhs_h_current end if energy%h(c) = energy%h_old(c) + params%dt * rhs_h_update energy_update_delta_rate_local = energy_update_delta_rate_local + & params%rho * (energy%h(c) - energy%h_old(c)) * mesh%cells(c)%volume / params%dt energy_update_rhs_integral_local = energy_update_rhs_integral_local + & params%rho * rhs_h_update * mesh%cells(c)%volume operator_consistent_rho_h_integral_local = operator_consistent_rho_h_integral_local + & params%rho * energy%h(c) * mesh%cells(c)%volume if (allocated(energy%operator_consistent_rho_h)) then energy%operator_consistent_rho_h(c) = params%rho * energy%h(c) end if end if energy%rhs_old(c) = rhs_h_current end do call profiler_stop('Energy_Flux_Update') energy%rhs_old_valid = .true. update_send(1) = energy_update_delta_rate_local update_send(2) = energy_update_rhs_integral_local update_send(3) = operator_consistent_rho_h_integral_local call MPI_Allreduce(update_send, update_recv, 3, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr) if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing direct energy update closure') energy_update_delta_rate_global = update_recv(1) energy_update_rhs_integral_global = update_recv(2) operator_consistent_rho_h_integral_global = update_recv(3) energy_update_balance_global = energy_update_delta_rate_global - energy_update_rhs_integral_global energy_update_denom = max(abs(energy_update_delta_rate_global), abs(energy_update_rhs_integral_global), tiny_safe) energy%last_energy_update_delta_rate_integral = energy_update_delta_rate_global energy%last_energy_update_rhs_integral = energy_update_rhs_integral_global energy%last_energy_update_balance_defect = energy_update_balance_global energy%relative_last_energy_update_balance_defect = abs(energy_update_balance_global) / energy_update_denom energy%cumulative_energy_update_delta_integral = & energy%cumulative_energy_update_delta_integral + params%dt * energy_update_delta_rate_global energy%cumulative_energy_update_rhs_integral = & energy%cumulative_energy_update_rhs_integral + params%dt * energy_update_rhs_integral_global energy%last_operator_consistent_rho_h_integral = operator_consistent_rho_h_integral_global energy%operator_consistent_rho_h_available = 1 energy%last_conductive_boundary_flux_out = conductive_boundary_flux_out_local energy%last_conductive_boundary_flux_available = 1 qrad_integral_step_local = zero rho_species_hdiff_integral_step_local = zero do c = flow%first_cell, flow%last_cell vol = mesh%cells(c)%volume if (allocated(energy%qrad)) qrad_integral_step_local = qrad_integral_step_local + energy%qrad(c) * vol if (allocated(energy%species_enthalpy_diffusion)) then if (use_variable_density_energy) then rho_diag = max(transport%rho(c), tiny_safe) else rho_diag = params%rho end if rho_species_hdiff_integral_step_local = rho_species_hdiff_integral_step_local + & rho_diag * energy%species_enthalpy_diffusion(c) * vol end if end do energy%cumulative_boundary_rho_h_advective_flux_out = & energy%cumulative_boundary_rho_h_advective_flux_out + params%dt * advective_boundary_rho_h_flux_out_local energy%cumulative_boundary_rho_h_conductive_flux_out = & energy%cumulative_boundary_rho_h_conductive_flux_out + params%dt * conductive_boundary_flux_out_local energy%cumulative_qrad_integral = energy%cumulative_qrad_integral + params%dt * qrad_integral_step_local energy%cumulative_rho_species_hdiff_integral = & energy%cumulative_rho_species_hdiff_integral + params%dt * rho_species_hdiff_integral_step_local energy%cumulative_energy_budget_available = 1 call write_variable_density_energy_debug(mesh, flow, params, fields, energy, & 'post_flux_pre_cantera_postsync', species_Y, transport, step) if (allocated(grad_h)) deallocate(grad_h) if (associated(hk_species)) deallocate(hk_species) if (associated(hk_boundary)) deallocate(hk_boundary) if (allocated(sp_diff_flux)) deallocate(sp_diff_flux) if (allocated(sp_Y_face_lin)) deallocate(sp_Y_face_lin) if (allocated(Y_point)) deallocate(Y_point) 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_owned(mesh, flow, 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. ! Pack related scalar fields to reduce neighbor-message latency. call profiler_start('Energy_Final_Exchange') if (allocated(energy%lambda) .and. allocated(energy%species_enthalpy_diffusion)) then call flow_exchange_cell_scalars(flow, energy%h, energy%T, energy%lambda, energy%species_enthalpy_diffusion) else if (allocated(energy%lambda)) then call flow_exchange_cell_scalars(flow, energy%h, energy%T, energy%lambda) else if (allocated(energy%species_enthalpy_diffusion)) then call flow_exchange_cell_scalars(flow, energy%h, energy%T, energy%species_enthalpy_diffusion) else call flow_exchange_cell_scalars(flow, energy%h, energy%T) end if call profiler_stop('Energy_Final_Exchange') end subroutine advance_energy_transport