Update the guarded variable-density low-Mach divergence source.
Constant-density mode returns without changing the zero source. In
variable-density mode the next projection targets
with the conservative source
The advective term is evaluated from corrected volumetric face fluxes as
a cell-local outward sum,
.
Interior values are centered; physical boundary faces use a
zero-normal-density-gradient assumption because no explicit density
boundary state exists yet. The routine updates fields%divergence_source
for the next projection and advances transport%rho_old after doing so.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(mesh_t), | intent(in) | :: | mesh |
Computational mesh context. |
||
| type(flow_mpi_t), | intent(inout) | :: | flow |
MPI parallelization communicator data. |
||
| type(bc_set_t), | intent(in) | :: | bc |
Boundary condition set. |
||
| type(case_params_t), | intent(in) | :: | params |
Global simulation configuration parameters. |
||
| type(transport_properties_t), | intent(inout) | :: | transport |
Thermodynamic property array containing density. |
||
| type(flow_fields_t), | intent(inout) | :: | fields |
Flow fields where the divergence source will be populated. |
||
| real(kind=rk), | intent(in), | optional | :: | species_Y(:,:) | ||
| real(kind=rk), | intent(in), | optional | :: | T_state(:) |
subroutine update_lowmach_divergence_source_from_density(mesh, flow, bc, params, transport, fields, species_Y, T_state) 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(transport_properties_t), intent(inout) :: transport type(flow_fields_t), intent(inout) :: fields real(rk), intent(in), optional :: species_Y(:,:) real(rk), intent(in), optional :: T_state(:) integer :: c !< Cell index counter. integer :: lf !< Local cell-to-face loop counter. integer :: fid !< Global face index. integer :: nb !< Neighbor cell index across the face. integer :: first_cell !< Starting cell index for the MPI partition. integer :: last_cell !< Ending cell index for the MPI partition. real(rk) :: vol !< Cell volume [m^3]. real(rk) :: rho_c !< Active cell density at current time level [kg/m^3]. real(rk) :: rho_old_c !< Cell density at the previous time level [kg/m^3]. real(rk) :: rho_nb !< Neighbor cell density [kg/m^3]. real(rk) :: rho_f !< Face density resolved by upwind/interpolation [kg/m^3]. real(rk) :: flux_out !< Outward volumetric face flux [m^3/s]. real(rk) :: history_source !< Temporal history source contribution (drho/dt) / rho. real(rk) :: advective_density_term !< Sum of outward convective density flux [kg/s]. real(rk) :: advective_source !< Advective source contribution -(u . grad(rho)) / rho. real(rk) :: rho_floor !< Small positive density safety floor to prevent division-by-zero. logical :: have_face_flux !< Flag indicating whether the face volumetric flux field is allocated. ! Local boundary-evaluation variables real(rk) :: T_b !< Interpolated boundary temperature [K]. real(rk) :: Y_b(max_species) !< Interpolated boundary species mass fraction vector. real(rk) :: h_sens_b !< Boundary enthalpy placeholder. real(rk) :: rho_b !< Resolved thermodynamic density at boundary [kg/m^3]. real(rk) :: cp_b !< Boundary specific heat placeholder. real(rk) :: lambda_b !< Boundary conductivity placeholder. real(rk) :: mu_b !< Boundary dynamic viscosity placeholder. logical :: is_dirichlet_T !< True if temperature boundary is Dirichlet. logical :: is_dirichlet_Y !< True if species boundaries are Dirichlet. integer :: k !< Species index loop counter. if (.not. params%enable_variable_density) return if (.not. allocated(transport%rho)) return if (.not. allocated(fields%divergence_source)) then allocate(fields%divergence_source(mesh%ncells)) fields%divergence_source = 0.0_rk else if (size(fields%divergence_source) /= mesh%ncells) then deallocate(fields%divergence_source) allocate(fields%divergence_source(mesh%ncells)) fields%divergence_source = 0.0_rk end if if (.not. allocated(transport%rho_old)) then allocate(transport%rho_old(mesh%ncells)) transport%rho_old = transport%rho else if (size(transport%rho_old) /= mesh%ncells) then deallocate(transport%rho_old) allocate(transport%rho_old(mesh%ncells)) transport%rho_old = transport%rho end if if (size(transport%rho) < mesh%ncells) return rho_floor = tiny(1.0_rk) have_face_flux = allocated(fields%face_flux) if (have_face_flux) have_face_flux = size(fields%face_flux) >= mesh%nfaces first_cell = max(1, flow%first_cell) last_cell = min(mesh%ncells, flow%last_cell) if (last_cell < first_cell) return do c = first_cell, last_cell vol = mesh%cells(c)%volume rho_c = max(transport%rho(c), rho_floor) rho_old_c = transport%rho_old(c) if (params%dt > rho_floor) then history_source = (rho_old_c - rho_c) / (rho_c * params%dt) else history_source = 0.0_rk end if advective_density_term = 0.0_rk if (have_face_flux .and. vol > rho_floor) then do lf = 1, mesh%ncell_faces(c) fid = mesh%cell_faces(lf, c) if (mesh%faces(fid)%owner == c) then flux_out = fields%face_flux(fid) nb = mesh%faces(fid)%neighbor if (nb <= 0) nb = mesh%faces(fid)%periodic_neighbor else flux_out = -fields%face_flux(fid) nb = mesh%faces(fid)%owner end if if (nb > 0 .and. nb <= mesh%ncells) then rho_nb = max(transport%rho(nb), rho_floor) if (params%upwind_dilatation_density) then if (flux_out >= 0.0_rk) then rho_f = rho_c else rho_f = rho_nb end if else rho_f = 0.5_rk * (rho_c + rho_nb) end if else ! Physical boundary face. If variable density and Cantera are enabled, ! and species/energy are present, evaluate the boundary density using Cantera. rho_f = rho_c if (params%enable_variable_density .and. present(species_Y) .and. present(T_state)) then if (allocated(transport%rho) .and. size(T_state) >= c) then ! 1. Retrieve the boundary temperature call boundary_temperature(mesh, bc, fid, T_state(c), T_b, is_dirichlet_T) ! 2. Retrieve the boundary species mass fractions Y_b = zero do k = 1, params%nspecies call boundary_species(mesh, bc, fid, k, species_Y(k, c), Y_b(k), is_dirichlet_Y) end do ! 3. Evaluate the boundary state density using Cantera call evaluate_cantera_boundary_state(params, T_b, Y_b(1:params%nspecies), & h_sens_b, rho_b, cp_b, lambda_b, mu_b) if (params%upwind_dilatation_density) then if (flux_out >= 0.0_rk) then rho_f = rho_c else rho_f = rho_b end if else rho_f = 0.5_rk * (rho_c + rho_b) end if end if end if end if advective_density_term = advective_density_term + flux_out * (rho_f - rho_c) end do advective_density_term = advective_density_term / vol end if advective_source = -advective_density_term / rho_c fields%divergence_source(c) = history_source + advective_source ! Advance density history only after storing the source that will be ! used by the next projection. transport%rho_old(c) = transport%rho(c) end do call flow_exchange_cell_scalar(flow, fields%divergence_source) call flow_exchange_cell_scalar(flow, transport%rho_old) end subroutine update_lowmach_divergence_source_from_density