write_species_energy_conservation_diagnostics

public subroutine write_species_energy_conservation_diagnostics(mesh, flow, params, fields, species, energy, transport, step, time)

Write aggregate species and enthalpy conservation trend diagnostics.
This diagnostics-only routine is a first-pass conservation tracker for the variable-density path. It reports global integrals, boundary-flux estimates, and step-to-step balance defects for:
total mass
transported species mass
transported species sum
rho*h
Species boundary fluxes use stored face mass flux and owner-cell species values on physical boundary faces. The rho*h flux is advective only; the aggregate enthalpy defect therefore excludes reconstructed conductive boundary fluxes and should be interpreted as a trend diagnostic, not a complete energy closure proof.
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(species_fields_t),	intent(in)	::	species
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_species_energy_conservation_diagnostics(mesh, flow, params, fields, species, 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(species_fields_t), intent(in) :: species
      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, f, k, owner, nb, ierr, unit_id, nsp
      logical, save :: aggregate_initialized = .false.
      logical, save :: species_initialized = .false.
      logical, save :: previous_initialized = .false.
      character(len=1024) :: filename
      character(len=64) :: label
      real(rk) :: vol, rho_c, y_sum, h_c, t_c, mflux
      real(rk) :: dt_since_previous, denom
      real(rk) :: local_volume, volume_total
      real(rk) :: local_mass, total_mass
      real(rk) :: local_boundary_mass_flux, net_boundary_mass_flux
      real(rk) :: local_species_mass_sum, species_mass_sum
      real(rk) :: local_boundary_species_flux_sum, net_boundary_species_flux_sum
      real(rk) :: local_y_sum_min, y_sum_min
      real(rk) :: local_y_sum_max, y_sum_max
      real(rk) :: local_y_sum_int, y_sum_int
      real(rk) :: local_y_sum_l2_num, y_sum_l2_num
      real(rk) :: species_sum_mean, species_sum_l2
      real(rk) :: local_y_sum_minus_one_abs_max, y_sum_minus_one_abs_max
      real(rk) :: local_y_sum_minus_one_l2_num, y_sum_minus_one_l2_num
      real(rk) :: species_sum_minus_one_l2
      real(rk) :: local_rho_h_integral, rho_h_integral
      real(rk) :: local_boundary_rho_h_advective_flux, net_boundary_rho_h_advective_flux
      real(rk) :: local_qrad_integral, qrad_integral
      real(rk) :: local_rho_species_hdiff_integral, rho_species_hdiff_integral
      real(rk) :: local_h_min, h_min
      real(rk) :: local_h_max, h_max
      real(rk) :: local_t_min, t_min
      real(rk) :: local_t_max, t_max
      real(rk) :: delta_mass, mass_rate, mass_balance_defect, rel_mass_balance_defect
      real(rk) :: delta_rho_h, rho_h_rate, rho_h_advective_balance_defect_no_conduction
      real(rk), save :: previous_time = 0.0_rk
      real(rk), save :: previous_mass = 0.0_rk
      real(rk), save :: previous_rho_h = 0.0_rk
      real(rk), allocatable, save :: previous_species_mass(:)
      real(rk), allocatable :: local_species_mass(:), species_mass(:)
      real(rk), allocatable :: local_boundary_species_flux(:), boundary_species_flux(:)
      real(rk), allocatable :: species_delta_mass(:), species_mass_rate(:)
      real(rk), allocatable :: species_balance_defect(:), species_relative_balance_defect(:)
      real(rk) :: send_val, recv_val

      if (.not. params%write_diagnostics) return

      nsp = 0
      if (params%enable_species) nsp = max(0, species%nspecies)

      allocate(local_species_mass(max(1, nsp)))
      allocate(species_mass(max(1, nsp)))
      allocate(local_boundary_species_flux(max(1, nsp)))
      allocate(boundary_species_flux(max(1, nsp)))
      allocate(species_delta_mass(max(1, nsp)))
      allocate(species_mass_rate(max(1, nsp)))
      allocate(species_balance_defect(max(1, nsp)))
      allocate(species_relative_balance_defect(max(1, nsp)))

      local_species_mass = 0.0_rk
      species_mass = 0.0_rk
      local_boundary_species_flux = 0.0_rk
      boundary_species_flux = 0.0_rk
      species_delta_mass = 0.0_rk
      species_mass_rate = 0.0_rk
      species_balance_defect = 0.0_rk
      species_relative_balance_defect = 0.0_rk

      local_volume = 0.0_rk
      local_mass = 0.0_rk
      local_boundary_mass_flux = 0.0_rk
      local_species_mass_sum = 0.0_rk
      local_boundary_species_flux_sum = 0.0_rk
      local_y_sum_min = huge(1.0_rk)
      local_y_sum_max = -huge(1.0_rk)
      local_y_sum_int = 0.0_rk
      local_y_sum_l2_num = 0.0_rk
      local_y_sum_minus_one_abs_max = 0.0_rk
      local_y_sum_minus_one_l2_num = 0.0_rk
      local_rho_h_integral = 0.0_rk
      local_boundary_rho_h_advective_flux = 0.0_rk
      local_qrad_integral = 0.0_rk
      local_rho_species_hdiff_integral = 0.0_rk
      local_h_min = huge(1.0_rk)
      local_h_max = -huge(1.0_rk)
      local_t_min = huge(1.0_rk)
      local_t_max = -huge(1.0_rk)

      do c = 1, mesh%ncells
         if (.not. flow%owned(c)) cycle

         vol = mesh%cells(c)%volume
         if (vol <= 0.0_rk) cycle

         rho_c = params%rho
         if (allocated(transport%rho)) then
            if (size(transport%rho) >= c) rho_c = transport%rho(c)
         end if

         local_volume = local_volume + vol
         local_mass = local_mass + rho_c * vol

         y_sum = 0.0_rk
         if (nsp > 0 .and. allocated(species%Y)) then
            if (size(species%Y, 2) >= c) then
               do k = 1, nsp
                  if (size(species%Y, 1) >= k) then
                     local_species_mass(k) = local_species_mass(k) + rho_c * species%Y(k, c) * vol
                     y_sum = y_sum + species%Y(k, c)
                  end if
               end do
            end if
         end if

         if (nsp > 0) then
            local_species_mass_sum = local_species_mass_sum + rho_c * y_sum * vol
            local_y_sum_min = min(local_y_sum_min, y_sum)
            local_y_sum_max = max(local_y_sum_max, y_sum)
            local_y_sum_int = local_y_sum_int + y_sum * vol
            local_y_sum_l2_num = local_y_sum_l2_num + y_sum * y_sum * vol
            local_y_sum_minus_one_abs_max = max(local_y_sum_minus_one_abs_max, abs(y_sum - 1.0_rk))
            local_y_sum_minus_one_l2_num = local_y_sum_minus_one_l2_num + (y_sum - 1.0_rk)**2 * vol
         end if

         if (params%enable_energy .and. allocated(energy%h)) then
            if (size(energy%h) >= c) then
               h_c = energy%h(c)
               local_rho_h_integral = local_rho_h_integral + rho_c * h_c * vol
               local_h_min = min(local_h_min, h_c)
               local_h_max = max(local_h_max, h_c)
            end if
         end if

         if (params%enable_energy .and. allocated(energy%T)) then
            if (size(energy%T) >= c) then
               t_c = energy%T(c)
               local_t_min = min(local_t_min, t_c)
               local_t_max = max(local_t_max, t_c)
            end if
         end if

         if (params%enable_energy .and. allocated(energy%qrad)) then
            if (size(energy%qrad) >= c) local_qrad_integral = local_qrad_integral + energy%qrad(c) * vol
         end if

         if (params%enable_energy .and. allocated(energy%species_enthalpy_diffusion)) then
            if (size(energy%species_enthalpy_diffusion) >= c) then
               local_rho_species_hdiff_integral = local_rho_species_hdiff_integral + &
                  rho_c * energy%species_enthalpy_diffusion(c) * vol
            end if
         end if
      end do

      if (allocated(fields%mass_flux)) then
         do f = 1, mesh%nfaces
            owner = mesh%faces(f)%owner
            nb = mesh%faces(f)%neighbor
            if (nb <= 0) nb = mesh%faces(f)%periodic_neighbor

            if (nb <= 0 .and. owner >= 1 .and. owner <= mesh%ncells) then
               if (.not. flow%owned(owner)) cycle

               mflux = fields%mass_flux(f)
               local_boundary_mass_flux = local_boundary_mass_flux + mflux

               if (nsp > 0 .and. allocated(species%Y)) then
                  if (size(species%Y, 2) >= owner) then
                     do k = 1, nsp
                        if (size(species%Y, 1) >= k) then
                           local_boundary_species_flux(k) = local_boundary_species_flux(k) + mflux * species%Y(k, owner)
                           local_boundary_species_flux_sum = local_boundary_species_flux_sum + mflux * species%Y(k, owner)
                        end if
                     end do
                  end if
               end if

               if (params%enable_energy .and. allocated(energy%h)) then
                  if (size(energy%h) >= owner) then
                     local_boundary_rho_h_advective_flux = local_boundary_rho_h_advective_flux + mflux * energy%h(owner)
                  end if
               end if
            end if
         end do
      end if

      call reduce_sum(local_volume, volume_total, flow)
      call reduce_sum(local_mass, total_mass, flow)
      call reduce_sum(local_boundary_mass_flux, net_boundary_mass_flux, flow)
      call reduce_sum(local_species_mass_sum, species_mass_sum, flow)
      call reduce_sum(local_boundary_species_flux_sum, net_boundary_species_flux_sum, flow)
      call reduce_sum(local_y_sum_int, y_sum_int, flow)
      call reduce_sum(local_y_sum_l2_num, y_sum_l2_num, flow)
      call reduce_sum(local_y_sum_minus_one_l2_num, y_sum_minus_one_l2_num, flow)
      call reduce_sum(local_rho_h_integral, rho_h_integral, flow)
      call reduce_sum(local_boundary_rho_h_advective_flux, net_boundary_rho_h_advective_flux, flow)
      call reduce_sum(local_qrad_integral, qrad_integral, flow)
      call reduce_sum(local_rho_species_hdiff_integral, rho_species_hdiff_integral, flow)

      if (nsp > 0) then
         send_val = local_y_sum_min
         call MPI_Allreduce(send_val, recv_val, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
         if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing species sum min')
         y_sum_min = recv_val

         send_val = local_y_sum_max
         call MPI_Allreduce(send_val, recv_val, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
         if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing species sum max')
         y_sum_max = recv_val

         send_val = local_y_sum_minus_one_abs_max
         call MPI_Allreduce(send_val, recv_val, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
         if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing species sum minus one max')
         y_sum_minus_one_abs_max = recv_val
      else
         y_sum_min = 0.0_rk
         y_sum_max = 0.0_rk
         y_sum_minus_one_abs_max = 0.0_rk
      end if

      if (params%enable_energy .and. allocated(energy%h)) then
         send_val = local_h_min
         call MPI_Allreduce(send_val, recv_val, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
         if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing h min')
         h_min = recv_val

         send_val = local_h_max
         call MPI_Allreduce(send_val, recv_val, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
         if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing h max')
         h_max = recv_val
      else
         h_min = 0.0_rk
         h_max = 0.0_rk
      end if

      if (params%enable_energy .and. allocated(energy%T)) then
         send_val = local_t_min
         call MPI_Allreduce(send_val, recv_val, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
         if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing T min')
         t_min = recv_val

         send_val = local_t_max
         call MPI_Allreduce(send_val, recv_val, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
         if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing T max')
         t_max = recv_val
      else
         t_min = 0.0_rk
         t_max = 0.0_rk
      end if

      if (nsp > 0) then
         call MPI_Allreduce(local_species_mass, species_mass, nsp, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
         if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing species mass integrals')

         call MPI_Allreduce(local_boundary_species_flux, boundary_species_flux, nsp, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
         if (ierr /= MPI_SUCCESS) call fatal_error('output', 'MPI failure reducing species boundary fluxes')
      end if

      if (volume_total > 0.0_rk) then
         species_sum_mean = y_sum_int / volume_total
         species_sum_l2 = sqrt(max(0.0_rk, y_sum_l2_num / volume_total))
         species_sum_minus_one_l2 = sqrt(max(0.0_rk, y_sum_minus_one_l2_num / volume_total))
      else
         species_sum_mean = 0.0_rk
         species_sum_l2 = 0.0_rk
         species_sum_minus_one_l2 = 0.0_rk
      end if

      if (previous_initialized) then
         dt_since_previous = time - previous_time
      else
         dt_since_previous = 0.0_rk
      end if

      if (previous_initialized .and. dt_since_previous > tiny(1.0_rk)) then
         delta_mass = total_mass - previous_mass
         mass_rate = delta_mass / dt_since_previous
         mass_balance_defect = mass_rate + net_boundary_mass_flux
         denom = max(abs(net_boundary_mass_flux), abs(mass_rate), tiny(1.0_rk))
         rel_mass_balance_defect = abs(mass_balance_defect) / denom

         delta_rho_h = rho_h_integral - previous_rho_h
         rho_h_rate = delta_rho_h / dt_since_previous
         rho_h_advective_balance_defect_no_conduction = rho_h_rate + net_boundary_rho_h_advective_flux - &
                                                        qrad_integral - rho_species_hdiff_integral

         if (nsp > 0) then
            if (.not. allocated(previous_species_mass)) allocate(previous_species_mass(nsp))
            if (size(previous_species_mass) /= nsp) then
               deallocate(previous_species_mass)
               allocate(previous_species_mass(nsp))
               previous_species_mass = species_mass
            end if

            do k = 1, nsp
               species_delta_mass(k) = species_mass(k) - previous_species_mass(k)
               species_mass_rate(k) = species_delta_mass(k) / dt_since_previous
               species_balance_defect(k) = species_mass_rate(k) + boundary_species_flux(k)
               denom = max(abs(species_mass_rate(k)), abs(boundary_species_flux(k)), tiny(1.0_rk))
               species_relative_balance_defect(k) = abs(species_balance_defect(k)) / denom
            end do
         end if
      else
         delta_mass = 0.0_rk
         mass_rate = 0.0_rk
         mass_balance_defect = 0.0_rk
         rel_mass_balance_defect = 0.0_rk
         delta_rho_h = 0.0_rk
         rho_h_rate = 0.0_rk
         rho_h_advective_balance_defect_no_conduction = 0.0_rk
      end if

      if (flow%rank == 0) then
         call execute_command_line('mkdir -p ' // trim(params%output_dir) // '/diagnostics')

         filename = trim(params%output_dir) // '/diagnostics/species_energy_conservation.csv'
         if (.not. aggregate_initialized) then
            open(newunit=unit_id, file=trim(filename), status='replace', action='write')
            aggregate_initialized = .true.
            write(unit_id,'(a)') 'step,time,volume_total,total_mass,net_boundary_mass_flux,' // &
               'delta_time_since_previous,delta_mass_since_previous,mass_rate_since_previous,' // &
               'mass_balance_defect_since_previous,relative_mass_balance_defect_since_previous,' // &
               'transported_species_mass_sum,net_boundary_species_mass_flux_sum,' // &
               'species_sum_min,species_sum_max,species_sum_mean,species_sum_l2,' // &
               'species_sum_minus_one_abs_max,species_sum_minus_one_l2,' // &
               'rho_h_integral,net_boundary_rho_h_advective_flux,' // &
               'qrad_integral,rho_species_enthalpy_diffusion_integral,' // &
               'delta_rho_h_since_previous,rho_h_rate_since_previous,' // &
               'rho_h_advective_balance_defect_no_conduction,h_min,h_max,T_min,T_max'
         else
            open(newunit=unit_id, file=trim(filename), status='unknown', position='append', action='write')
         end if

         call write_csv_value_i(unit_id, step)
         call write_csv_value_r(unit_id, time)
         call write_csv_value_r(unit_id, volume_total)
         call write_csv_value_r(unit_id, total_mass)
         call write_csv_value_r(unit_id, net_boundary_mass_flux)
         call write_csv_value_r(unit_id, dt_since_previous)
         call write_csv_value_r(unit_id, delta_mass)
         call write_csv_value_r(unit_id, mass_rate)
         call write_csv_value_r(unit_id, mass_balance_defect)
         call write_csv_value_r(unit_id, rel_mass_balance_defect)
         call write_csv_value_r(unit_id, species_mass_sum)
         call write_csv_value_r(unit_id, net_boundary_species_flux_sum)
         call write_csv_value_r(unit_id, y_sum_min)
         call write_csv_value_r(unit_id, y_sum_max)
         call write_csv_value_r(unit_id, species_sum_mean)
         call write_csv_value_r(unit_id, species_sum_l2)
         call write_csv_value_r(unit_id, y_sum_minus_one_abs_max)
         call write_csv_value_r(unit_id, species_sum_minus_one_l2)
         call write_csv_value_r(unit_id, rho_h_integral)
         call write_csv_value_r(unit_id, net_boundary_rho_h_advective_flux)
         call write_csv_value_r(unit_id, qrad_integral)
         call write_csv_value_r(unit_id, rho_species_hdiff_integral)
         call write_csv_value_r(unit_id, delta_rho_h)
         call write_csv_value_r(unit_id, rho_h_rate)
         call write_csv_value_r(unit_id, rho_h_advective_balance_defect_no_conduction)
         call write_csv_value_r(unit_id, h_min)
         call write_csv_value_r(unit_id, h_max)
         call write_csv_value_r(unit_id, t_min)
         call write_csv_value_r(unit_id, t_max)
         write(unit_id,*)
         close(unit_id)

         filename = trim(params%output_dir) // '/diagnostics/species_integrals.csv'
         if (nsp > 0) then
            if (.not. species_initialized) then
               open(newunit=unit_id, file=trim(filename), status='replace', action='write')
               species_initialized = .true.
               write(unit_id,'(a)', advance='no') 'step,time'
               do k = 1, nsp
                  write(label,'("Y",i0)') k
                  call write_csv_header(unit_id, trim(label)//'_mass_integral')
                  call write_csv_header(unit_id, trim(label)//'_net_boundary_mass_flux')
                  call write_csv_header(unit_id, trim(label)//'_delta_mass_since_previous')
                  call write_csv_header(unit_id, trim(label)//'_mass_rate_since_previous')
                  call write_csv_header(unit_id, trim(label)//'_balance_defect_since_previous')
                  call write_csv_header(unit_id, trim(label)//'_relative_balance_defect_since_previous')
               end do
               write(unit_id,*)
            else
               open(newunit=unit_id, file=trim(filename), status='unknown', position='append', action='write')
            end if

            call write_csv_value_i(unit_id, step)
            call write_csv_value_r(unit_id, time)
            do k = 1, nsp
               call write_csv_value_r(unit_id, species_mass(k))
               call write_csv_value_r(unit_id, boundary_species_flux(k))
               call write_csv_value_r(unit_id, species_delta_mass(k))
               call write_csv_value_r(unit_id, species_mass_rate(k))
               call write_csv_value_r(unit_id, species_balance_defect(k))
               call write_csv_value_r(unit_id, species_relative_balance_defect(k))
            end do
            write(unit_id,*)
            close(unit_id)
         end if
      end if

      previous_time = time
      previous_mass = total_mass
      previous_rho_h = rho_h_integral
      previous_initialized = .true.

      if (nsp > 0) then
         if (.not. allocated(previous_species_mass)) allocate(previous_species_mass(nsp))
         if (size(previous_species_mass) /= nsp) then
            deallocate(previous_species_mass)
            allocate(previous_species_mass(nsp))
         end if
         previous_species_mass = species_mass(1:nsp)
      end if

      deallocate(local_species_mass)
      deallocate(species_mass)
      deallocate(local_boundary_species_flux)
      deallocate(boundary_species_flux)
      deallocate(species_delta_mass)
      deallocate(species_mass_rate)
      deallocate(species_balance_defect)
      deallocate(species_relative_balance_defect)

   end subroutine write_species_energy_conservation_diagnostics
write_species_energy_conservation_diagnostics Subroutine

Contents

Source Code

public subroutine write_species_energy_conservation_diagnostics(mesh, flow, params, fields, species, energy, transport, step, time)

Arguments

Source Code
