write_variable_density_boundary_residual_scan

private subroutine write_variable_density_boundary_residual_scan(mesh, flow, params, fields, energy, transport, step, time)

Scan all owned boundary-adjacent cells for low-Mach projection residuals.
A boundary-adjacent cell is any owned cell touching at least one physical boundary face, identified by neighbor <= 0 and periodic_neighbor <= 0. The scan recomputes div(u) and div(rho u) from the current face fluxes and writes per-rank boundary-cell rows plus a rank-0 global summary.
Arguments

Type	Intent		Name
type(mesh_t),	intent(in)	::	mesh
type(flow_mpi_t),	intent(in)	::	flow
type(case_params_t),	intent(in)	::	params
type(flow_fields_t),	intent(in)	::	fields
type(energy_fields_t),	intent(in)	::	energy
type(transport_properties_t),	intent(in)	::	transport
integer,	intent(in)	::	step
real(kind=rk),	intent(in)	::	time
Source Code

   subroutine write_variable_density_boundary_residual_scan(mesh, flow, params, fields, energy, transport, step, time)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(case_params_t), intent(in) :: params
      type(flow_fields_t), intent(in) :: fields
      type(energy_fields_t), intent(in) :: energy
      type(transport_properties_t), intent(in) :: transport
      integer, intent(in) :: step
      real(rk), intent(in) :: time

      integer :: c, lf, fid, nb, ierr, unit_id
      integer :: boundary_face_count
      integer :: local_boundary_count, local_interior_count
      integer :: boundary_count, interior_count
      integer :: local_boundary_worst_cell, local_interior_worst_cell
      integer :: local_boundary_winner_rank, global_boundary_winner_rank
      integer :: local_interior_winner_rank, global_interior_winner_rank
      integer :: global_boundary_worst_cell, global_interior_worst_cell
      integer :: mpi_int_send(1), mpi_int_recv(1)
      logical :: is_boundary_cell
      logical, save :: vd_boundary_scan_initialized = .false.
      character(len=32) :: rank_suffix
      character(len=1024) :: filename
      real(rk) :: vol, rho_c, s_c, div_c, div_res, mdot_div_c
      real(rk) :: flux_out, mdot_out, mass_flux_value
      real(rk) :: T_c, h_c
      real(rk) :: local_boundary_volume, local_interior_volume
      real(rk) :: boundary_volume, interior_volume
      real(rk) :: local_boundary_res_max, local_boundary_res_l2_num
      real(rk) :: local_interior_res_max, local_interior_res_l2_num
      real(rk) :: boundary_res_max, boundary_res_l2_num, boundary_res_l2
      real(rk) :: interior_res_max, interior_res_l2_num, interior_res_l2
      real(rk) :: local_boundary_div_min, local_boundary_div_max
      real(rk) :: local_boundary_S_min, local_boundary_S_max
      real(rk) :: local_boundary_mdot_min, local_boundary_mdot_max
      real(rk) :: boundary_div_min, boundary_div_max
      real(rk) :: boundary_S_min, boundary_S_max
      real(rk) :: boundary_mdot_min, boundary_mdot_max
      real(rk) :: local_boundary_worst_divu, local_boundary_worst_S
      real(rk) :: local_boundary_worst_mdot_div, local_boundary_worst_rho
      real(rk) :: local_boundary_worst_T, local_boundary_worst_h
      real(rk) :: local_interior_worst_divu, local_interior_worst_S
      real(rk) :: local_interior_worst_mdot_div, local_interior_worst_rho
      real(rk) :: local_interior_worst_T, local_interior_worst_h
      real(rk) :: global_boundary_worst_divu, global_boundary_worst_S
      real(rk) :: global_boundary_worst_mdot_div, global_boundary_worst_rho
      real(rk) :: global_boundary_worst_T, global_boundary_worst_h
      real(rk) :: global_interior_worst_divu, global_interior_worst_S
      real(rk) :: global_interior_worst_mdot_div, global_interior_worst_rho
      real(rk) :: global_interior_worst_T, global_interior_worst_h
      real(rk) :: mpi_reduce_send(1), mpi_reduce_recv(1)

      if (.not. params%enable_variable_density) return
      if (.not. params%write_diagnostics) return
      if (.not. allocated(transport%rho)) return

      write(rank_suffix, '(i0)') flow%rank

      filename = trim(params%output_dir) // '/diagnostics/variable_density_boundary_residual_cells_rank' // &
                 trim(adjustl(rank_suffix)) // '.csv'
      if (.not. vd_boundary_scan_initialized) then
         open(newunit=unit_id, file=trim(filename), status='replace', action='write')
         write(unit_id,'(a)') 'step,time,rank,cell,boundary_face_count,volume,rho,T,h,S,divu,divu_minus_S,mass_flux_div,abs_residual'
      else
         open(newunit=unit_id, file=trim(filename), status='unknown', position='append', action='write')
      end if

      local_boundary_count = 0
      local_interior_count = 0
      local_boundary_volume = 0.0_rk
      local_interior_volume = 0.0_rk
      local_boundary_res_max = 0.0_rk
      local_interior_res_max = 0.0_rk
      local_boundary_res_l2_num = 0.0_rk
      local_interior_res_l2_num = 0.0_rk

      local_boundary_div_min = huge(1.0_rk)
      local_boundary_div_max = -huge(1.0_rk)
      local_boundary_S_min = huge(1.0_rk)
      local_boundary_S_max = -huge(1.0_rk)
      local_boundary_mdot_min = huge(1.0_rk)
      local_boundary_mdot_max = -huge(1.0_rk)

      local_boundary_worst_cell = -1
      local_interior_worst_cell = -1
      local_boundary_worst_divu = 0.0_rk
      local_boundary_worst_S = 0.0_rk
      local_boundary_worst_mdot_div = 0.0_rk
      local_boundary_worst_rho = 0.0_rk
      local_boundary_worst_T = 0.0_rk
      local_boundary_worst_h = 0.0_rk
      local_interior_worst_divu = 0.0_rk
      local_interior_worst_S = 0.0_rk
      local_interior_worst_mdot_div = 0.0_rk
      local_interior_worst_rho = 0.0_rk
      local_interior_worst_T = 0.0_rk
      local_interior_worst_h = 0.0_rk

      do c = flow%first_cell, flow%last_cell
         vol = mesh%cells(c)%volume
         rho_c = transport%rho(c)

         s_c = 0.0_rk
         if (allocated(fields%divergence_source)) s_c = fields%divergence_source(c)

         T_c = 0.0_rk
         h_c = 0.0_rk
         if (allocated(energy%T)) T_c = energy%T(c)
         if (allocated(energy%h)) h_c = energy%h(c)

         div_c = 0.0_rk
         mdot_div_c = 0.0_rk
         boundary_face_count = 0

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

            nb = mesh%faces(fid)%neighbor
            if (nb <= 0) then
               if (mesh%faces(fid)%periodic_neighbor <= 0) boundary_face_count = boundary_face_count + 1
            end if

            if (allocated(fields%face_flux)) then
               if (mesh%faces(fid)%owner == c) then
                  flux_out = fields%face_flux(fid)
               else
                  flux_out = -fields%face_flux(fid)
               end if
               div_c = div_c + flux_out / vol
            end if

            if (allocated(fields%mass_flux)) then
               mass_flux_value = fields%mass_flux(fid)
               if (mesh%faces(fid)%owner == c) then
                  mdot_out = mass_flux_value
               else
                  mdot_out = -mass_flux_value
               end if
               mdot_div_c = mdot_div_c + mdot_out / vol
            end if
         end do

         is_boundary_cell = boundary_face_count > 0
         div_res = div_c - s_c

         if (is_boundary_cell) then
            local_boundary_count = local_boundary_count + 1
            local_boundary_volume = local_boundary_volume + vol
            local_boundary_res_l2_num = local_boundary_res_l2_num + div_res * div_res * vol
            local_boundary_div_min = min(local_boundary_div_min, div_c)
            local_boundary_div_max = max(local_boundary_div_max, div_c)
            local_boundary_S_min = min(local_boundary_S_min, s_c)
            local_boundary_S_max = max(local_boundary_S_max, s_c)
            local_boundary_mdot_min = min(local_boundary_mdot_min, mdot_div_c)
            local_boundary_mdot_max = max(local_boundary_mdot_max, mdot_div_c)

            write(unit_id,'(i0,",",ES26.16E4,",",i0,",",i0,",",i0,9(",",ES26.16E4))') &
                 step, time, flow%rank, c, boundary_face_count, vol, rho_c, T_c, h_c, s_c, &
                 div_c, div_res, mdot_div_c, abs(div_res)

            if (abs(div_res) > local_boundary_res_max) then
               local_boundary_res_max = abs(div_res)
               local_boundary_worst_cell = c
               local_boundary_worst_divu = div_c
               local_boundary_worst_S = s_c
               local_boundary_worst_mdot_div = mdot_div_c
               local_boundary_worst_rho = rho_c
               local_boundary_worst_T = T_c
               local_boundary_worst_h = h_c
            end if
         else
            local_interior_count = local_interior_count + 1
            local_interior_volume = local_interior_volume + vol
            local_interior_res_l2_num = local_interior_res_l2_num + div_res * div_res * vol

            if (abs(div_res) > local_interior_res_max) then
               local_interior_res_max = abs(div_res)
               local_interior_worst_cell = c
               local_interior_worst_divu = div_c
               local_interior_worst_S = s_c
               local_interior_worst_mdot_div = mdot_div_c
               local_interior_worst_rho = rho_c
               local_interior_worst_T = T_c
               local_interior_worst_h = h_c
            end if
         end if
      end do

      close(unit_id)

      if (local_boundary_count == 0) then
         local_boundary_div_min = 0.0_rk
         local_boundary_div_max = 0.0_rk
         local_boundary_S_min = 0.0_rk
         local_boundary_S_max = 0.0_rk
         local_boundary_mdot_min = 0.0_rk
         local_boundary_mdot_max = 0.0_rk
      end if

      mpi_int_send(1) = local_boundary_count
      call MPI_Allreduce(mpi_int_send, mpi_int_recv, 1, MPI_INTEGER, MPI_SUM, flow%comm, ierr)
      boundary_count = mpi_int_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing boundary cell count')

      mpi_int_send(1) = local_interior_count
      call MPI_Allreduce(mpi_int_send, mpi_int_recv, 1, MPI_INTEGER, MPI_SUM, flow%comm, ierr)
      interior_count = mpi_int_recv(1)
      if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing interior cell count')

      mpi_reduce_send(1) = local_boundary_volume
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      boundary_volume = mpi_reduce_recv(1)

      mpi_reduce_send(1) = local_interior_volume
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      interior_volume = mpi_reduce_recv(1)

      mpi_reduce_send(1) = local_boundary_res_max
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      boundary_res_max = mpi_reduce_recv(1)

      mpi_reduce_send(1) = local_interior_res_max
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      interior_res_max = mpi_reduce_recv(1)

      mpi_reduce_send(1) = local_boundary_res_l2_num
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      boundary_res_l2_num = mpi_reduce_recv(1)

      mpi_reduce_send(1) = local_interior_res_l2_num
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      interior_res_l2_num = mpi_reduce_recv(1)

      mpi_reduce_send(1) = local_boundary_div_min
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
      boundary_div_min = mpi_reduce_recv(1)

      mpi_reduce_send(1) = local_boundary_div_max
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      boundary_div_max = mpi_reduce_recv(1)

      mpi_reduce_send(1) = local_boundary_S_min
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
      boundary_S_min = mpi_reduce_recv(1)

      mpi_reduce_send(1) = local_boundary_S_max
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      boundary_S_max = mpi_reduce_recv(1)

      mpi_reduce_send(1) = local_boundary_mdot_min
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
      boundary_mdot_min = mpi_reduce_recv(1)

      mpi_reduce_send(1) = local_boundary_mdot_max
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      boundary_mdot_max = mpi_reduce_recv(1)

      if (boundary_volume > 0.0_rk) then
         boundary_res_l2 = sqrt(boundary_res_l2_num / boundary_volume)
      else
         boundary_res_l2 = 0.0_rk
      end if

      if (interior_volume > 0.0_rk) then
         interior_res_l2 = sqrt(interior_res_l2_num / interior_volume)
      else
         interior_res_l2 = 0.0_rk
      end if

      if (abs(local_boundary_res_max - boundary_res_max) <= max(1.0e-12_rk, 1.0e-10_rk * max(1.0_rk, boundary_res_max))) then
         local_boundary_winner_rank = flow%rank
      else
         local_boundary_winner_rank = huge(1)
      end if

      mpi_int_send(1) = local_boundary_winner_rank
      call MPI_Allreduce(mpi_int_send, mpi_int_recv, 1, MPI_INTEGER, MPI_MIN, flow%comm, ierr)
      global_boundary_winner_rank = mpi_int_recv(1)

      if (abs(local_interior_res_max - interior_res_max) <= max(1.0e-12_rk, 1.0e-10_rk * max(1.0_rk, interior_res_max))) then
         local_interior_winner_rank = flow%rank
      else
         local_interior_winner_rank = huge(1)
      end if

      mpi_int_send(1) = local_interior_winner_rank
      call MPI_Allreduce(mpi_int_send, mpi_int_recv, 1, MPI_INTEGER, MPI_MIN, flow%comm, ierr)
      global_interior_winner_rank = mpi_int_recv(1)

      if (flow%rank == global_boundary_winner_rank) then
         mpi_int_send(1) = local_boundary_worst_cell
      else
         mpi_int_send(1) = 0
      end if
      call MPI_Allreduce(mpi_int_send, mpi_int_recv, 1, MPI_INTEGER, MPI_SUM, flow%comm, ierr)
      global_boundary_worst_cell = mpi_int_recv(1)

      if (flow%rank == global_interior_winner_rank) then
         mpi_int_send(1) = local_interior_worst_cell
      else
         mpi_int_send(1) = 0
      end if
      call MPI_Allreduce(mpi_int_send, mpi_int_recv, 1, MPI_INTEGER, MPI_SUM, flow%comm, ierr)
      global_interior_worst_cell = mpi_int_recv(1)

      mpi_reduce_send(1) = merge(local_boundary_worst_divu, 0.0_rk, flow%rank == global_boundary_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_boundary_worst_divu = mpi_reduce_recv(1)

      mpi_reduce_send(1) = merge(local_boundary_worst_S, 0.0_rk, flow%rank == global_boundary_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_boundary_worst_S = mpi_reduce_recv(1)

      mpi_reduce_send(1) = merge(local_boundary_worst_mdot_div, 0.0_rk, flow%rank == global_boundary_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_boundary_worst_mdot_div = mpi_reduce_recv(1)

      mpi_reduce_send(1) = merge(local_boundary_worst_rho, 0.0_rk, flow%rank == global_boundary_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_boundary_worst_rho = mpi_reduce_recv(1)

      mpi_reduce_send(1) = merge(local_boundary_worst_T, 0.0_rk, flow%rank == global_boundary_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_boundary_worst_T = mpi_reduce_recv(1)

      mpi_reduce_send(1) = merge(local_boundary_worst_h, 0.0_rk, flow%rank == global_boundary_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_boundary_worst_h = mpi_reduce_recv(1)

      mpi_reduce_send(1) = merge(local_interior_worst_divu, 0.0_rk, flow%rank == global_interior_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_interior_worst_divu = mpi_reduce_recv(1)

      mpi_reduce_send(1) = merge(local_interior_worst_S, 0.0_rk, flow%rank == global_interior_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_interior_worst_S = mpi_reduce_recv(1)

      mpi_reduce_send(1) = merge(local_interior_worst_mdot_div, 0.0_rk, flow%rank == global_interior_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_interior_worst_mdot_div = mpi_reduce_recv(1)

      mpi_reduce_send(1) = merge(local_interior_worst_rho, 0.0_rk, flow%rank == global_interior_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_interior_worst_rho = mpi_reduce_recv(1)

      mpi_reduce_send(1) = merge(local_interior_worst_T, 0.0_rk, flow%rank == global_interior_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_interior_worst_T = mpi_reduce_recv(1)

      mpi_reduce_send(1) = merge(local_interior_worst_h, 0.0_rk, flow%rank == global_interior_winner_rank)
      call MPI_Allreduce(mpi_reduce_send, mpi_reduce_recv, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      global_interior_worst_h = mpi_reduce_recv(1)

      if (flow%rank == 0) then
         filename = trim(params%output_dir) // '/diagnostics/variable_density_boundary_residual_summary.csv'
         if (.not. vd_boundary_scan_initialized) then
            open(newunit=unit_id, file=trim(filename), status='replace', action='write')
            write(unit_id,'(a)') 'step,time,boundary_cell_count,interior_cell_count,boundary_res_max,boundary_res_l2,interior_res_max,interior_res_l2,boundary_divu_min,boundary_divu_max,boundary_S_min,boundary_S_max,boundary_mass_flux_div_min,boundary_mass_flux_div_max,boundary_worst_rank,boundary_worst_cell,boundary_worst_divu,boundary_worst_S,boundary_worst_mass_flux_div,boundary_worst_rho,boundary_worst_T,boundary_worst_h,interior_worst_rank,interior_worst_cell,interior_worst_divu,interior_worst_S,interior_worst_mass_flux_div,interior_worst_rho,interior_worst_T,interior_worst_h'
         else
            open(newunit=unit_id, file=trim(filename), status='unknown', position='append', action='write')
         end if

         write(unit_id,'(i0,",",ES26.16E4,2(",",i0),10(",",ES26.16E4),2(",",i0),6(",",ES26.16E4),2(",",i0),6(",",ES26.16E4))') &
              step, time, boundary_count, interior_count, boundary_res_max, boundary_res_l2, &
              interior_res_max, interior_res_l2, boundary_div_min, boundary_div_max, boundary_S_min, boundary_S_max, &
              boundary_mdot_min, boundary_mdot_max, global_boundary_winner_rank, global_boundary_worst_cell, &
              global_boundary_worst_divu, global_boundary_worst_S, global_boundary_worst_mdot_div, &
              global_boundary_worst_rho, global_boundary_worst_T, global_boundary_worst_h, &
              global_interior_winner_rank, global_interior_worst_cell, global_interior_worst_divu, &
              global_interior_worst_S, global_interior_worst_mdot_div, global_interior_worst_rho, &
              global_interior_worst_T, global_interior_worst_h
         close(unit_id)
      end if

      vd_boundary_scan_initialized = .true.

   end subroutine write_variable_density_boundary_residual_scan
write_variable_density_boundary_residual_scan Subroutine

Contents

Source Code

private subroutine write_variable_density_boundary_residual_scan(mesh, flow, params, fields, energy, transport, step, time)

Arguments

Source Code
