write_projection_audit_one_cell Subroutine

private subroutine write_projection_audit_one_cell(mesh, flow, fields, energy, transport, step, time, c, audit_label, unit_cells, unit_faces)

Write one audited cell and all its faces.

Arguments

Type IntentOptional Attributes Name
type(mesh_t), intent(in) :: mesh
type(flow_mpi_t), intent(in) :: flow
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
integer, intent(in) :: c
character(len=*), intent(in) :: audit_label
integer, intent(in) :: unit_cells
integer, intent(in) :: unit_faces

Source Code

   subroutine write_projection_audit_one_cell(mesh, flow, fields, energy, transport, step, time, &
                                              c, audit_label, unit_cells, unit_faces)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      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, intent(in) :: c
      character(len=*), intent(in) :: audit_label
      integer, intent(in) :: unit_cells
      integer, intent(in) :: unit_faces

      integer :: lf, fid, nb, periodic_nb, owner_cell
      integer :: boundary_face_count, physical_boundary_face_count
      integer :: is_boundary_face, is_physical_boundary_face
      real(rk) :: vol, rho_c, T_c, h_c, s_c, div_c, div_res, mdot_div_c
      real(rk) :: required_flux_sum, actual_flux_sum, flux_defect, mass_flux_sum
      real(rk) :: cell_pressure, owner_pressure, neighbor_pressure, dp_neighbor_minus_cell
      real(rk) :: face_flux_value, outward_face_flux, mass_flux_value, outward_mass_flux

      if (c < 1) return

      vol = mesh%cells(c)%volume
      rho_c = transport%rho(c)
      T_c = 0.0_rk
      h_c = 0.0_rk
      s_c = 0.0_rk

      if (allocated(energy%T)) T_c = energy%T(c)
      if (allocated(energy%h)) h_c = energy%h(c)
      if (allocated(fields%divergence_source)) s_c = fields%divergence_source(c)

      call compute_projection_audit_cell_balance(mesh, fields, c, div_c, mdot_div_c, &
                                                 boundary_face_count, physical_boundary_face_count)

      div_res = div_c - s_c
      required_flux_sum = s_c * vol
      actual_flux_sum = div_c * vol
      flux_defect = actual_flux_sum - required_flux_sum
      mass_flux_sum = mdot_div_c * vol
      cell_pressure = 0.0_rk
      if (allocated(fields%p)) then
         if (c >= lbound(fields%p, 1) .and. c <= ubound(fields%p, 1)) cell_pressure = fields%p(c)
      end if

      write(unit_cells,'(i0,",",ES26.16E4,",",i0,",",a,",",i0,12(",",ES26.16E4),2(",",i0),",",ES26.16E4)') &
           step, time, flow%rank, trim(audit_label), c, vol, rho_c, T_c, h_c, s_c, div_c, div_res, &
           mdot_div_c, required_flux_sum, actual_flux_sum, flux_defect, mass_flux_sum, &
           boundary_face_count, physical_boundary_face_count, cell_pressure

      do lf = 1, mesh%ncell_faces(c)
         fid = mesh%cell_faces(lf, c)
         nb = mesh%faces(fid)%neighbor
         periodic_nb = mesh%faces(fid)%periodic_neighbor

         is_boundary_face = 0
         is_physical_boundary_face = 0
         if (nb <= 0) then
            is_boundary_face = 1
            if (periodic_nb <= 0) is_physical_boundary_face = 1
         end if

         face_flux_value = 0.0_rk
         if (allocated(fields%face_flux)) face_flux_value = fields%face_flux(fid)

         mass_flux_value = 0.0_rk
         if (allocated(fields%mass_flux)) mass_flux_value = fields%mass_flux(fid)

         if (mesh%faces(fid)%owner == c) then
            outward_face_flux = face_flux_value
            outward_mass_flux = mass_flux_value
         else
            outward_face_flux = -face_flux_value
            outward_mass_flux = -mass_flux_value
         end if

         owner_cell = mesh%faces(fid)%owner
         owner_pressure = 0.0_rk
         neighbor_pressure = 0.0_rk
         if (allocated(fields%p)) then
            if (owner_cell >= lbound(fields%p, 1) .and. owner_cell <= ubound(fields%p, 1)) then
               owner_pressure = fields%p(owner_cell)
            end if
            if (nb >= lbound(fields%p, 1) .and. nb <= ubound(fields%p, 1)) then
               neighbor_pressure = fields%p(nb)
            end if
         end if
         dp_neighbor_minus_cell = neighbor_pressure - cell_pressure

         write(unit_faces,'(i0,",",ES26.16E4,",",i0,",",a,8(",",i0),8(",",ES26.16E4))') &
              step, time, flow%rank, trim(audit_label), c, lf, fid, mesh%faces(fid)%owner, nb, periodic_nb, &
              is_boundary_face, is_physical_boundary_face, owner_pressure, cell_pressure, neighbor_pressure, &
              dp_neighbor_minus_cell, face_flux_value, outward_face_flux, mass_flux_value, outward_mass_flux
      end do

   end subroutine write_projection_audit_one_cell