balance_neumann_outlet_flux Subroutine

private subroutine balance_neumann_outlet_flux(mesh, flow, bc, face_flux)

Adjusts flux at Neumann outlets to ensure strict global mass balance.

This is critical for solvers without a pressure-pinned cell, as numerical drift in the Poisson RHS can lead to a singular system.

Arguments

Type IntentOptional Attributes Name
type(mesh_t), intent(in) :: mesh
type(flow_mpi_t), intent(in) :: flow
type(bc_set_t), intent(in) :: bc
real(kind=rk), intent(inout) :: face_flux(:)

Calls

proc~~balance_neumann_outlet_flux~~CallsGraph proc~balance_neumann_outlet_flux mod_flow_projection::balance_neumann_outlet_flux mpi_allreduce mpi_allreduce proc~balance_neumann_outlet_flux->mpi_allreduce proc~check_mpi mod_flow_projection::check_mpi proc~balance_neumann_outlet_flux->proc~check_mpi proc~fatal_error mod_kinds::fatal_error proc~balance_neumann_outlet_flux->proc~fatal_error proc~patch_type_for_face mod_bc::patch_type_for_face proc~balance_neumann_outlet_flux->proc~patch_type_for_face proc~profiler_start mod_profiler::profiler_start proc~balance_neumann_outlet_flux->proc~profiler_start proc~profiler_stop mod_profiler::profiler_stop proc~balance_neumann_outlet_flux->proc~profiler_stop proc~check_mpi->proc~fatal_error 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

Called by

proc~~balance_neumann_outlet_flux~~CalledByGraph proc~balance_neumann_outlet_flux mod_flow_projection::balance_neumann_outlet_flux proc~advance_projection_step mod_flow_projection::advance_projection_step proc~advance_projection_step->proc~balance_neumann_outlet_flux program~lowmach_react_hex lowmach_react_hex program~lowmach_react_hex->proc~advance_projection_step

Source Code

   subroutine balance_neumann_outlet_flux(mesh, flow, bc, face_flux)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(bc_set_t), intent(in) :: bc
      real(rk), intent(inout) :: face_flux(:)

      integer :: f, owner, btype, ierr
      real(rk) :: local_net_flux
      real(rk) :: global_net_flux
      real(rk) :: local_outlet_area
      real(rk) :: global_outlet_area
      real(rk) :: correction_per_area

      if (pressure_cache%has_dirichlet_pressure) return

      ! If the case has no Neumann outlet faces, there is nothing to balance.
      ! This avoids two unnecessary MPI_Allreduce calls per timestep for closed
      ! domains such as lid-driven cavity.
      if (.not. pressure_cache%has_neumann_outlet) return

      local_net_flux = zero
      local_outlet_area = zero

      do f = 1, mesh%nfaces
         if (mesh%faces(f)%neighbor /= 0) cycle

         btype = patch_type_for_face(mesh, bc, f)

         if (btype == bc_periodic) cycle

         owner = mesh%faces(f)%owner
         if (.not. flow%owned(owner)) cycle

         local_net_flux = local_net_flux + face_flux(f)

         if (btype == bc_neumann) then
            local_outlet_area = local_outlet_area + mesh%faces(f)%area
         end if
      end do

      call profiler_start('MPI_Communication')
      call MPI_Allreduce(local_net_flux, global_net_flux, 1, &
                         MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      call check_mpi(ierr, 'outlet balance net flux')

      call MPI_Allreduce(local_outlet_area, global_outlet_area, 1, &
                         MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      call check_mpi(ierr, 'outlet balance outlet area')
      call profiler_stop('MPI_Communication')

      if (global_outlet_area <= tiny_safe) then
         if (abs(global_net_flux) > tiny_safe) then
            call fatal_error('flow', 'nonzero boundary flux but no neumann outlet faces found')
         end if
         return
      end if

      correction_per_area = -global_net_flux / global_outlet_area

      do f = 1, mesh%nfaces
         if (mesh%faces(f)%neighbor /= 0) cycle

         btype = patch_type_for_face(mesh, bc, f)
         if (btype /= bc_neumann) cycle
         owner = mesh%faces(f)%owner
         if (.not. flow%owned(owner)) cycle

         face_flux(f) = face_flux(f) + correction_per_area * mesh%faces(f)%area
      end do
   end subroutine balance_neumann_outlet_flux