lowmach_react_hex Program

program lowmach_react_hex
   use mpi_f08
   use mod_kinds, only : rk, zero, one, output_unit, fatal_error
   use mod_input, only : case_params_t, read_case_params
   use mod_mesh_types, only : mesh_t, mesh_finalize
   use mod_mesh_io, only : read_native_mesh
   use mod_mpi_flow, only : flow_mpi_t, mpi_flow_startup, mpi_flow_shutdown, &
                            flow_mpi_initialize, flow_mpi_finalize
   use mod_mpi_radiation, only : radiation_mpi_t, radiation_mpi_initialize, &
                                 radiation_mpi_finalize
   use mod_bc, only : bc_set_t, build_bc_set, finalize_bc_set
   use mod_fields, only : flow_fields_t, initialize_fields, finalize_fields
   use mod_flow_projection, only : solver_stats_t, advance_projection_step, compute_flow_diagnostics, compute_and_update_cfl
   use mod_output, only : prepare_output, write_diagnostics_header, write_diagnostics_row, &
                          write_vtu_unstructured, write_mesh_summary, write_pvd_collection
   use mod_species, only : species_fields_t, initialize_species, finalize_species, &
                           advance_species_transport
   use mod_energy, only : energy_fields_t, initialize_energy, finalize_energy, &
                           advance_energy_transport, &
                           write_energy_diagnostics_header, write_energy_diagnostics_row
   use mod_transport_properties, only : transport_properties_t, initialize_transport, &
                                        finalize_transport, update_transport_properties
   use mod_profiler, only : profiler_start, profiler_stop, profiler_report, profiler_configure
   implicit none

   type(case_params_t) :: params    !< Parsed simulation parameters.
   type(mesh_t) :: mesh             !< Global computational mesh.
   type(flow_mpi_t) :: flow_mpi     !< MPI context for flow solver.
   type(radiation_mpi_t) :: rad_mpi !< MPI context for radiation solver.
   type(bc_set_t) :: bc             !< Boundary condition definitions.
   type(flow_fields_t) :: fields    !< Hydrodynamic fields (U, P).
   type(solver_stats_t) :: stats    !< Runtime diagnostics and solver performance metrics.
   type(species_fields_t) :: species !< Species mass fractions.
   type(energy_fields_t) :: energy   !< Enthalpy/temperature/radiation-source fields.
   type(transport_properties_t) :: transport !< Physical properties (rho, mu, Dk).

   character(len=256) :: case_file
   integer :: step, ierr, world_rank
   logical :: do_output_step
   real(rk) :: time, t0

   ! 1. Initialize MPI and read runtime configuration.
   call mpi_flow_startup()
   call MPI_Comm_rank(MPI_COMM_WORLD, world_rank, ierr)
   if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI_Comm_rank failed')

   t0 = real(mpi_wtime(), rk)

   call get_case_filename(case_file)
   call read_case_params(trim(case_file), params)
   call profiler_configure(params%enable_profiling, params%nested_profiling)

   if (world_rank == 0) then
      write(output_unit,'(a,l1,a,l1)') 'profiling config: enable=', &
         params%enable_profiling, ' nested=', params%nested_profiling
   end if

   call profiler_start('Total_Simulation')

   ! 2. Setup simulation environment and read input data.
   call read_native_mesh(trim(params%mesh_dir), mesh)
   
   ! 3. Initialize parallel contexts.
   call initialize_transport(mesh, params, transport)
   call flow_mpi_initialize(mesh, flow_mpi, MPI_COMM_WORLD, max(4, params%nspecies))
   call radiation_mpi_initialize(rad_mpi, MPI_COMM_WORLD)
   
   ! 4. Prepare physical fields and BCs.
   call build_bc_set(mesh, params, bc)
   call initialize_fields(mesh, fields)
   if (params%enable_species) then
      call initialize_species(mesh, params, species)
   end if
   if (params%enable_energy) then
      if (params%enable_species) then
         call initialize_energy(mesh, params, energy, species%Y)
      else
         call initialize_energy(mesh, params, energy)
      end if
   end if

   if (flow_mpi%rank == 0 .and. params%enable_species) then
      print *, "species: ", species%nspecies
   end if

   ! 5. Setup output files.
   call prepare_output(params, flow_mpi)
   call write_mesh_summary(params, flow_mpi, mesh)
   call write_diagnostics_header(params, flow_mpi)
   call write_energy_diagnostics_header(params, flow_mpi)
   call write_pvd_collection(params, flow_mpi, params%nsteps, params%output_interval, params%dt)

   ! 6. Initial diagnostics and visualization snapshot.
   time = zero
   call compute_flow_diagnostics(mesh, flow_mpi, bc, params, fields, stats)
   call write_diagnostics_row(params, flow_mpi, 0, time, stats)
   call write_energy_diagnostics_row(mesh, flow_mpi, params, energy, 0, time)
   call write_vtu_unstructured(params, flow_mpi, mesh, fields, species, energy, transport, 0)

   if (flow_mpi%rank == 0) then
      write(output_unit,'(a)') 'LowMachReact-Hex hexahedral low-Mach FV solver'
      write(output_unit,'(a,a)') 'case: ', trim(case_file)
      write(output_unit,'(a,i0)') 'cells: ', mesh%ncells
      write(output_unit,'(a,i0)') 'flow/chemistry MPI ranks: ', flow_mpi%nprocs
      write(output_unit,'(a,i0)') 'radiation MPI ranks: ', rad_mpi%nprocs
      write(output_unit,'(a,es12.5,a)') 'Flow density mode: constant rho = ', params%rho, ' kg/m^3'
      if (params%enable_variable_nu) then
         write(output_unit,'(a)') 'Flow viscosity mode: variable Cantera mu/nu enabled'
      else
         write(output_unit,'(a,es12.5,a)') 'Flow viscosity mode: constant nu = ', params%nu, ' m^2/s'
      end if
      if (params%enable_energy .and. params%enable_cantera_thermo) then
         write(output_unit,'(a)') 'Cantera rho_thermo: diagnostic only, not used by projection'
      end if
      if (params%enable_cantera_fluid .or. params%enable_cantera_species) then
         write(output_unit,'(a,i0,a)') 'Cantera transport update interval: ', &
            params%transport_update_interval, ' step(s) [mu/D_k only]'
         if (params%enable_energy) then
            write(output_unit,'(a)') 'Cantera transport temperature source: energy%T'
         else
            write(output_unit,'(a,es12.5,a)') 'Cantera transport temperature source: background_temp = ', &
               params%background_temp, ' K'
         end if
      end if
      if (params%enable_energy .and. params%enable_cantera_thermo) then
         write(output_unit,'(a)') 'Cantera thermo update interval: every energy step'
         write(output_unit,'(a,es12.5,a)') 'Cantera thermodynamic pressure p0: ', &
            params%background_press, ' Pa'
      end if
   end if

   ! 7. Main Time-Stepping Loop
   do step = 1, params%nsteps
      do_output_step = (mod(step, params%output_interval) == 0) .or. step == params%nsteps

      ! A. Adaptive time-stepping / CFL diagnostics.
      ! Dynamic dt requires CFL every step. For fixed-dt runs, compute CFL only
      ! when diagnostics are written.
      if (params%use_dynamic_dt .or. do_output_step) then
         call profiler_start('CFL_Update')
         call compute_and_update_cfl(mesh, flow_mpi, params, fields, stats)
         call profiler_stop('CFL_Update')
      end if

      ! B. Dynamic transport-property update via Cantera.
      ! transport_update_interval gates Cantera mu/D_k refresh; enable_variable_nu controls whether mu affects flow; flow density remains params%rho.
      ! Cantera h/T thermo recovery remains inside the energy update path.
      if (mod(step-1, params%transport_update_interval) == 0 .or. step == 1) then
         call profiler_start('Transport_Update')
         if (params%enable_species) then
            if (params%enable_energy) then
               call update_transport_properties(mesh, flow_mpi, params, species%Y, transport, T_state=energy%T)
            else
               call update_transport_properties(mesh, flow_mpi, params, species%Y, transport)
            end if
         else
            if (params%enable_energy) then
               call update_transport_properties(mesh, flow_mpi, params, transport=transport, T_state=energy%T)
            else
               call update_transport_properties(mesh, flow_mpi, params, transport=transport)
            end if
         end if
         call profiler_stop('Transport_Update')
      end if
      
      ! C. Advance Momentum & Pressure (Projection Method).
      call profiler_start('Projection_Step')
      call advance_projection_step(mesh, flow_mpi, bc, params, transport, fields, stats)
      call profiler_stop('Projection_Step')
      
      ! D. Advance Species Transport.
      if (params%enable_species) then
         call profiler_start('Species_Transport')
         call advance_species_transport(mesh, flow_mpi, bc, params, fields, species, transport)
         call profiler_stop('Species_Transport')
      end if
      
      ! E. Advance passive sensible-enthalpy transport.
      if (params%enable_energy) then
         call profiler_start('Energy_Transport')
         if (params%enable_species) then
            call advance_energy_transport(mesh, flow_mpi, bc, params, fields, energy, species%Y)
         else
            call advance_energy_transport(mesh, flow_mpi, bc, params, fields, energy)
         end if
         call profiler_stop('Energy_Transport')
      end if
      
      time = time + params%dt

      ! E. Periodic Output and Diagnostics.
      if (do_output_step) then
         stats%wall_time = real(mpi_wtime(), rk) - t0
         call profiler_start('Flow_Diagnostics')
         call compute_global_observables(mesh, flow_mpi, fields, species, transport, stats)
         call profiler_stop('Flow_Diagnostics')

         call profiler_start('Diagnostics_Write_Flow')
         call write_diagnostics_row(params, flow_mpi, step, time, stats)
         call profiler_stop('Diagnostics_Write_Flow')

         if (params%enable_energy) then
            call profiler_start('Diagnostics_Write_Energy')
            call write_energy_diagnostics_row(mesh, flow_mpi, params, energy, step, time)
            call profiler_stop('Diagnostics_Write_Energy')
         end if

         call profiler_start('Output_Write_VTU')
         call write_vtu_unstructured(params, flow_mpi, mesh, fields, species, energy, transport, step)
         call profiler_stop('Output_Write_VTU')

         if (flow_mpi%rank == 0) then
            write(output_unit,'(a,i0,2x,a,es12.5,2x,a,f8.2,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,i0,2x,a,es12.5)') &
               'step=', step, 'time=', time, 'clock=', stats%wall_time, 'dt=', params%dt, 'max_div=', stats%max_divergence, &
               'cfl=', stats%cfl, 'ke=', stats%kinetic_energy, 'piter=', stats%pressure_iterations, '|U|max=', stats%max_velocity
         end if
      end if
   end do

   ! 8. Clean up and finalize simulation.
   call finalize_bc_set(bc)
   call finalize_fields(fields)
   if (params%enable_species) then
      call finalize_species(species)
   end if
   if (params%enable_energy) then
      call finalize_energy(energy)
   end if
   call finalize_transport(transport)
   call mesh_finalize(mesh)

   call profiler_stop('Total_Simulation')
   call profiler_report(flow_mpi%comm, flow_mpi%rank, flow_mpi%nprocs)

   call radiation_mpi_finalize(rad_mpi)
   call flow_mpi_finalize(flow_mpi)
   call mpi_flow_shutdown()

contains

   !> Computes global integral quantities for diagnostic reporting.
   !!
   !! Performs Allreduce operations to find the maximum velocity magnitude, 
   !! total system mass, and minimum species mass fraction.
   subroutine compute_global_observables(mesh, flow, fields, species, transport, stats)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(flow_fields_t), intent(in) :: fields
      type(species_fields_t), intent(in) :: species
      type(transport_properties_t), intent(in) :: transport
      type(solver_stats_t), intent(inout) :: stats

      integer :: c
      real(rk) :: local_max_vel, local_min_y, local_total_mass
      real(rk) :: global_max_vel, global_min_y, global_total_mass
      real(rk) :: u, v, w, vel_mag

      local_max_vel = zero
      local_min_y = huge(one)
      local_total_mass = zero

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

         u = fields%u(1, c)
         v = fields%u(2, c)
         w = fields%u(3, c)
         vel_mag = sqrt(u**2 + v**2 + w**2)
         if (vel_mag > local_max_vel) local_max_vel = vel_mag

         local_total_mass = local_total_mass + (transport%rho(c) * mesh%cells(c)%volume)

         if (species%nspecies > 0) then
            local_min_y = min(local_min_y, minval(species%Y(:, c)))
         end if
      end do

      call MPI_Allreduce(local_max_vel, global_max_vel, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm)
      call MPI_Allreduce(local_total_mass, global_total_mass, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm)
      
      if (species%nspecies > 0) then
         call MPI_Allreduce(local_min_y, global_min_y, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm)
      else
         global_min_y = zero
      end if

      stats%max_velocity = global_max_vel
      stats%total_mass = global_total_mass
      stats%min_species_y = global_min_y

   end subroutine compute_global_observables


   !> Parses command line arguments to find the case configuration file.
   !!
   !! If no argument is provided, defaults to `case.nml`.
   subroutine get_case_filename(filename)
      character(len=*), intent(out) :: filename

      integer :: argc

      argc = command_argument_count()
      if (argc >= 1) then
         call get_command_argument(1, filename)
      else
         filename = 'case.nml'
      end if

      if (len_trim(filename) == 0) call fatal_error('main', 'empty case filename')
   end subroutine get_case_filename

end program lowmach_react_hex