m_global_parameters.fpp.f90

Go to the documentation of this file.

# 1 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
!>
!! @file
!! @brief Contains module m_global_parameters
 
# 1 "/home/runner/work/MFC/MFC/src/common/include/case.fpp" 1
! This file exists so that Fypp can be run without generating case.fpp files for
! each target. This is useful when generating documentation, for example. This
! should also let MFC be built with CMake directly, without invoking mfc.sh.
 
! For pre-process.
# 9 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
 
! For moving immersed boundaries in simulation
# 14 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
# 6 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp" 2
 
!> @brief Defines global parameters for the computational domain, simulation algorithm, and initial conditions
module m_global_parameters
 
#ifdef MFC_MPI
    use mpi  ! Message passing interface (MPI) module
#endif
 
    use m_derived_types  ! Definitions of the derived types
    use m_helper_basic  ! Functions to compare floating point numbers
    use m_thermochem, only: num_species
 
    implicit none
 
    ! Logistics
    integer                 :: num_procs                 !< Number of processors
    character(LEN=path_len) :: case_dir                  !< Case folder location
    logical                 :: old_grid                  !< Use existing grid data
    logical                 :: old_ic, non_axis_sym      !< Use existing IC data
    integer                 :: t_step_old, t_step_start  !< Existing IC/grid folder
    logical                 :: cfl_adap_dt, cfl_const_dt, cfl_dt
    integer                 :: n_start, n_start_old
 
    ! Computational Domain Parameters
 
    integer :: proc_rank  !< Rank of the local processor Number of cells in the x-, y- and z-coordinate directions
    integer :: m
    integer :: n
    integer :: p
 
    !> @name Max and min number of cells in a direction of each combination of x-,y-, and z-
    type(cell_num_bounds) :: cells_bounds
    integer(kind=8)       :: nglobal              !< Global number of cells in the domain
    integer               :: m_glb, n_glb, p_glb  !< Global number of cells in each direction
    integer               :: num_dims             !< Number of spatial dimensions
    integer               :: num_vels             !< Number of velocity components (different from num_dims for mhd)
    logical               :: cyl_coord
    integer               :: grid_geometry        !< Cylindrical coordinates (either axisymmetric or full 3D)
    !> Locations of cell-centers (cc) in x-, y- and z-directions, respectively
    real(wp), allocatable, dimension(:) :: x_cc, y_cc, z_cc
    !> Locations of cell-boundaries (cb) in x-, y- and z-directions, respectively
    real(wp), allocatable, dimension(:) :: x_cb, y_cb, z_cb
    real(wp) :: dx, dy, dz                             !< Minimum cell-widths in the x-, y- and z-coordinate directions
    type(bounds_info) :: x_domain, y_domain, z_domain  !< Locations of the domain bounds in the x-, y- and z-coordinate directions
    logical :: stretch_x, stretch_y, stretch_z         !< Grid stretching flags for the x-, y- and z-coordinate directions
    ! Grid stretching: a_x/a_y/a_z = rate, x_a/y_a/z_a = location
    real(wp) :: a_x, a_y, a_z
    integer  :: loops_x, loops_y, loops_z
    real(wp) :: x_a, y_a, z_a
    real(wp) :: x_b, y_b, z_b
 
    ! Simulation Algorithm Parameters
    integer            :: model_eqns                   !< Multicomponent flow model
    logical            :: relax                        !< activate phase change
    integer            :: relax_model                  !< Relax Model
    real(wp)           :: palpha_eps                   !< trigger parameter for the p relaxation procedure, phase change model
    real(wp)           :: ptgalpha_eps                 !< trigger parameter for the pTg relaxation procedure, phase change model
    integer            :: num_fluids                   !< Number of different fluids present in the flow
    logical            :: mpp_lim                      !< Alpha limiter
    integer            :: sys_size                     !< Number of unknowns in the system of equations
    integer            :: recon_type                   !< Reconstruction Type
    integer            :: weno_polyn                   !< Degree of the WENO polynomials (polyn)
    integer            :: muscl_polyn                  !< Degree of the MUSCL polynomials (polyn)
    integer            :: weno_order                   !< Order of accuracy for the WENO reconstruction
    integer            :: muscl_order                  !< Order of accuracy for the MUSCL reconstruction
    logical            :: hypoelasticity               !< activate hypoelasticity
    logical            :: hyperelasticity              !< activate hyperelasticity
    logical            :: elasticity                   !< elasticity modeling, true for hyper or hypo
    logical            :: mhd                          !< Magnetohydrodynamics
    logical            :: relativity                   !< Relativity for RMHD
    integer            :: b_size                       !< Number of components in the b tensor
    integer            :: tensor_size                  !< Number of components in the nonsymmetric tensor
    logical            :: pre_stress                   !< activate pre_stressed domain
    logical            :: cont_damage                  !< continuum damage modeling
    logical            :: hyper_cleaning               !< Hyperbolic cleaning for MHD
    logical            :: igr                          !< Use information geometric regularization
    integer            :: igr_order                    !< IGR reconstruction order
    logical, parameter :: chemistry = .false.  !< Chemistry modeling
    ! Annotations of the structure, i.e. the organization, of the state vectors
    type(eqn_idx_info)  :: eqn_idx   !< All conserved-variable equation index ranges and scalars.
    type(qbmm_idx_info) :: qbmm_idx  !< QBMM moment index mappings.
    ! Cell Indices for the (local) interior points (O-m, O-n, 0-p). Stands for "InDices With BUFFer".
    type(int_bounds_info) :: idwint(1:3)
 
    ! Cell indices (InDices With BUFFer): includes buffer except in pre_process
    type(int_bounds_info)      :: idwbuff(1:3)
    type(int_bounds_info)      :: bc_x, bc_y, bc_z       !< Boundary conditions in the x-, y- and z-coordinate directions
    integer                    :: shear_num              !< Number of shear stress components
    integer, dimension(3)      :: shear_indices          !< Indices of the stress components that represent shear stress
    integer                    :: shear_bc_flip_num      !< Number of shear stress components to reflect for boundary conditions
    integer, dimension(3, 2)   :: shear_bc_flip_indices  !< Shear stress BC reflection indices (1:3, 1:shear_BC_flip_num)
    logical                    :: parallel_io            !< Format of the data files
    logical                    :: file_per_process       !< type of data output
    integer                    :: precision              !< Precision of output files
    logical                    :: down_sample            !< Down-sample the output data
    logical                    :: mixlayer_vel_profile   !< Set hyperbolic tangent streamwise velocity profile
    real(wp)                   :: mixlayer_vel_coef      !< Coefficient for the hyperbolic tangent streamwise velocity profile
    logical                    :: mixlayer_perturb       !< Superimpose instability waves to surrounding fluid flow
    integer                    :: mixlayer_perturb_nk    !< Number of Fourier modes for perturbation with mixlayer_perturb flag
    real(wp)                   :: mixlayer_perturb_k0    !< Peak wavenumber for mixlayer perturbation (default: most unstable mode)
    logical                    :: simplex_perturb
    type(simplex_noise_params) :: simplex_params
    real(wp)                   :: pi_fac                 !< Factor for artificial pi_inf
    logical                    :: viscous
    logical                    :: bubbles_lagrange
 
    ! Perturb density of surrounding air so as to break symmetry of grid
    logical                             :: perturb_flow
    integer                             :: perturb_flow_fluid  !< Fluid to be perturbed with perturb_flow flag
    real(wp)                            :: perturb_flow_mag    !< Magnitude of perturbation with perturb_flow flag
    logical                             :: perturb_sph
    integer                             :: perturb_sph_fluid   !< Fluid to be perturbed with perturb_sph flag
    real(wp), dimension(num_fluids_max) :: fluid_rho
    logical                             :: elliptic_smoothing
    integer                             :: elliptic_smoothing_iters
    integer, allocatable, dimension(:)  :: proc_coords         !< Processor coordinates in MPI_CART_COMM
    integer, allocatable, dimension(:)  :: start_idx           !< Starting cell-center index of local processor in global grid
#ifdef MFC_MPI
    type(mpi_io_var), public :: mpi_io_data
    character(LEN=name_len)  :: mpiiofs
    integer                  :: mpi_info_int  !< MPI info for parallel IO with Lustre file systems
#endif
 
    ! Initial Condition Parameters
    integer                                                  :: num_patches     !< Number of patches composing initial condition
    type(ic_patch_parameters), dimension(num_patches_max)    :: patch_icpp      !< IC patch parameters (max: num_patches_max)
    integer                                                  :: num_bc_patches  !< Number of boundary condition patches
    logical                                                  :: bc_io           !< whether or not to save BC data
    type(bc_patch_parameters), dimension(num_bc_patches_max) :: patch_bc        !< Boundary condition patch parameters
 
    ! Fluids Physical Parameters
    type(physical_parameters), dimension(num_fluids_max) :: fluid_pp  !< Stiffened gas EOS parameters and Reynolds numbers per fluid
    ! Subgrid Bubble Parameters
    type(subgrid_bubble_physical_parameters) :: bub_pp
    real(wp)                                 :: rhoref, pref  !< Reference parameters for Tait EOS
    type(chemistry_parameters)               :: chem_params
    !> @name Bubble modeling
    !> @{
    integer                             :: nb
    real(wp)                            :: ca, web, re_inv, eu
    real(wp), dimension(:), allocatable :: weight, r0
    logical                             :: bubbles_euler
    logical                             :: qbmm  !< Quadrature moment method
    integer                             :: nmom  !< Number of carried moments
    real(wp)                            :: sigr, sigv, rhorv  !< standard deviations in R/V
    logical                             :: adv_n  !< Solve the number density equation and compute alpha from number density
    !> @}
 
    !> @name Immersed Boundaries
    !> @{
    logical                                                  :: ib        !< Turn immersed boundaries on
    integer                                                  :: num_ibs   !< Number of immersed boundaries
    integer                                                  :: np
    type(ib_patch_parameters), dimension(num_ib_patches_max) :: patch_ib  !< Immersed boundary patch parameters
    type(vec3_dt), allocatable, dimension(:)                 :: airfoil_grid_u, airfoil_grid_l
    !> @}
 
    !> @name Non-polytropic bubble gas compression
    !> @{
    logical :: polytropic
    logical :: polydisperse
    real(wp) :: poly_sigma
    integer :: dist_type  !< 1 = binormal, 2 = lognormal-normal
    integer :: thermal    !< 1 = adiabatic, 2 = isotherm, 3 = transfer
    real(wp) :: phi_vg, phi_gv, pe_c, tw, k_vl, k_gl
    real(wp) :: gam_m
    real(wp), dimension(:), allocatable :: pb0, mass_g0, mass_v0, pe_t, k_v, k_g
    real(wp), dimension(:), allocatable :: re_trans_t, re_trans_c, im_trans_t, im_trans_c, omegan
    real(wp) :: r0ref, p0ref, rho0ref, t0ref, ss, pv, vd, mu_l, mu_v, mu_g, gam_v, gam_g, m_v, m_g, cp_v, cp_g, r_v, r_g
    !> @}
 
    !> @name Surface Tension Modeling
    !> @{
    real(wp) :: sigma
    logical  :: surface_tension
    !> @}
 
    integer, allocatable, dimension(:,:,:) :: logic_grid
    type(pres_field)                       :: pb
    type(pres_field)                       :: mv
    real(wp)                               :: bx0        !< Constant magnetic field in the x-direction (1D)
    integer                                :: buff_size  !< Number of ghost cells for boundary condition storage
    logical                                :: fft_wrt
 
contains
 
    !> Assigns default values to user inputs prior to reading them in. This allows for an easier consistency check of these
    !! parameters once they are read from the input file.
    impure subroutine s_assign_default_values_to_user_inputs
 
        integer :: i  !< Generic loop operator
        ! Logistics
 
        case_dir = '.'
        old_grid = .false.
        old_ic = .false.
        t_step_old = dflt_int
        t_step_start = dflt_int
 
        cfl_adap_dt = .false.
        cfl_const_dt = .false.
        cfl_dt = .false.
        n_start = dflt_int
 
        ! Computational domain parameters
        m = dflt_int; n = 0; p = 0
 
        call s_update_cell_bounds(cells_bounds, m, n, p)
 
        cyl_coord = .false.
 
        x_domain%beg = dflt_real
        x_domain%end = dflt_real
        y_domain%beg = dflt_real
        y_domain%end = dflt_real
        z_domain%beg = dflt_real
        z_domain%end = dflt_real
 
        stretch_x = .false.
        stretch_y = .false.
        stretch_z = .false.
 
        a_x = dflt_real
        a_y = dflt_real
        a_z = dflt_real
        loops_x = 1
        loops_y = 1
        loops_z = 1
        x_a = dflt_real
        x_b = dflt_real
        y_a = dflt_real
        y_b = dflt_real
        z_a = dflt_real
        z_b = dflt_real
 
        ! Simulation algorithm parameters
        model_eqns = dflt_int
        relax = .false.
        relax_model = dflt_int
        palpha_eps = dflt_real
        ptgalpha_eps = dflt_real
        num_fluids = dflt_int
        recon_type = weno_type
        weno_order = dflt_int
        igr = .false.
        igr_order = dflt_int
        muscl_order = dflt_int
 
        hypoelasticity = .false.
        hyperelasticity = .false.
        elasticity = .false.
        pre_stress = .false.
        b_size = dflt_int
        tensor_size = dflt_int
        cont_damage = .false.
        hyper_cleaning = .false.
 
        mhd = .false.
        relativity = .false.
 
        bc_x%beg = dflt_int; bc_x%end = dflt_int
        bc_y%beg = dflt_int; bc_y%end = dflt_int
        bc_z%beg = dflt_int; bc_z%end = dflt_int
 
# 271 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
# 272 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
                bc_x%vb1 = 0._wp
                bc_x%ve1 = 0._wp
# 272 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
                bc_x%vb2 = 0._wp
                bc_x%ve2 = 0._wp
# 272 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
                bc_x%vb3 = 0._wp
                bc_x%ve3 = 0._wp
# 275 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
# 271 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
# 272 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
                bc_y%vb1 = 0._wp
                bc_y%ve1 = 0._wp
# 272 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
                bc_y%vb2 = 0._wp
                bc_y%ve2 = 0._wp
# 272 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
                bc_y%vb3 = 0._wp
                bc_y%ve3 = 0._wp
# 275 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
# 271 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
# 272 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
                bc_z%vb1 = 0._wp
                bc_z%ve1 = 0._wp
# 272 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
                bc_z%vb2 = 0._wp
                bc_z%ve2 = 0._wp
# 272 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
                bc_z%vb3 = 0._wp
                bc_z%ve3 = 0._wp
# 275 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
# 276 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
 
# 278 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
            bc_z%isothermal_in = .false.
            bc_z%isothermal_out = .false.
            bc_z%Twall_in = dflt_real
            bc_z%Twall_out = dflt_real
# 278 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
            bc_x%isothermal_in = .false.
            bc_x%isothermal_out = .false.
            bc_x%Twall_in = dflt_real
            bc_x%Twall_out = dflt_real
# 278 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
            bc_y%isothermal_in = .false.
            bc_y%isothermal_out = .false.
            bc_y%Twall_in = dflt_real
            bc_y%Twall_out = dflt_real
# 283 "/home/runner/work/MFC/MFC/src/pre_process/m_global_parameters.fpp"
 
        parallel_io = .false.
        file_per_process = .false.
        precision = 2
        down_sample = .false.
        viscous = .false.
        bubbles_lagrange = .false.
        mixlayer_vel_profile = .false.
        mixlayer_vel_coef = 1._wp
        mixlayer_perturb = .false.
        mixlayer_perturb_nk = 100
        mixlayer_perturb_k0 = 0.4446_wp
        perturb_flow = .false.
        perturb_flow_fluid = dflt_int
        perturb_flow_mag = dflt_real
        perturb_sph = .false.
        perturb_sph_fluid = dflt_int
        fluid_rho = dflt_real
        elliptic_smoothing_iters = dflt_int
        elliptic_smoothing = .false.
 
        fft_wrt = .false.
 
        simplex_perturb = .false.
        simplex_params%perturb_vel(:) = .false.
        simplex_params%perturb_vel_freq(:) = dflt_real
        simplex_params%perturb_vel_scale(:) = dflt_real
        simplex_params%perturb_vel_offset(:,:) = dflt_real
        simplex_params%perturb_dens(:) = .false.
        simplex_params%perturb_dens_freq(:) = dflt_real
        simplex_params%perturb_dens_scale(:) = dflt_real
        simplex_params%perturb_dens_offset(:,:) = dflt_real
 
        ! Initial condition parameters
        num_patches = dflt_int
 
        do i = 1, num_patches_max
            patch_icpp(i)%geometry = dflt_int
            patch_icpp(i)%model_scale(:) = 1._wp
            patch_icpp(i)%model_translate(:) = 0._wp
            patch_icpp(i)%model_filepath(:) = dflt_char
            patch_icpp(i)%model_spc = num_ray
            patch_icpp(i)%model_threshold = ray_tracing_threshold
            patch_icpp(i)%x_centroid = dflt_real
            patch_icpp(i)%y_centroid = dflt_real
            patch_icpp(i)%z_centroid = dflt_real
            patch_icpp(i)%length_x = dflt_real
            patch_icpp(i)%length_y = dflt_real
            patch_icpp(i)%length_z = dflt_real
            patch_icpp(i)%radius = dflt_real
            patch_icpp(i)%epsilon = dflt_real
            patch_icpp(i)%beta = dflt_real
            patch_icpp(i)%normal = dflt_real
            patch_icpp(i)%radii = dflt_real
            patch_icpp(i)%alter_patch = .false.
            patch_icpp(i)%alter_patch(0) = .true.
            patch_icpp(i)%smoothen = .false.
            patch_icpp(i)%smooth_patch_id = i
            patch_icpp(i)%smooth_coeff = dflt_real
            patch_icpp(i)%alpha_rho = dflt_real
            patch_icpp(i)%rho = dflt_real
            patch_icpp(i)%vel = dflt_real
            patch_icpp(i)%pres = dflt_real
            patch_icpp(i)%alpha = dflt_real
            patch_icpp(i)%gamma = dflt_real
            patch_icpp(i)%pi_inf = dflt_real
            patch_icpp(i)%cv = 0._wp
            patch_icpp(i)%qv = 0._wp
            patch_icpp(i)%qvp = 0._wp
            patch_icpp(i)%tau_e = 0._wp
            patch_icpp(i)%Bx = dflt_real
            patch_icpp(i)%By = dflt_real
            patch_icpp(i)%Bz = dflt_real
            patch_icpp(i)%a(2) = dflt_real
            patch_icpp(i)%a(3) = dflt_real
            patch_icpp(i)%a(4) = dflt_real
            patch_icpp(i)%a(5) = dflt_real
            patch_icpp(i)%a(6) = dflt_real
            patch_icpp(i)%a(7) = dflt_real
            patch_icpp(i)%a(8) = dflt_real
            patch_icpp(i)%a(9) = dflt_real
            patch_icpp(i)%non_axis_sym = .false.
            patch_icpp(i)%fourier_cos(:) = 0._wp
            patch_icpp(i)%fourier_sin(:) = 0._wp
            patch_icpp(i)%modal_clip_r_to_min = .false.
            patch_icpp(i)%modal_r_min = 1.e-12_wp
            patch_icpp(i)%modal_use_exp_form = .false.
            patch_icpp(i)%sph_har_coeff(:,:) = 0._wp
 
            ! should get all of r0's and v0's
            patch_icpp(i)%r0 = dflt_real
            patch_icpp(i)%v0 = dflt_real
 
            patch_icpp(i)%p0 = dflt_real
            patch_icpp(i)%m0 = dflt_real
 
            patch_icpp(i)%hcid = dflt_int
 
            if (chemistry) then
                patch_icpp(i)%Y(:) = 0._wp
            end if
        end do
 
        num_bc_patches = 0
        bc_io = .false.
 
        do i = 1, num_bc_patches_max
            patch_bc(i)%geometry = dflt_int
            patch_bc(i)%type = dflt_int
            patch_bc(i)%dir = dflt_int
            patch_bc(i)%loc = dflt_int
            patch_bc(i)%centroid(:) = dflt_real
            patch_bc(i)%length(:) = dflt_real
            patch_bc(i)%radius = dflt_real
        end do
 
        ! Tait EOS
        rhoref = dflt_real
        pref = dflt_real
 
        ! Bubble modeling
        bubbles_euler = .false.
        polytropic = .true.
        polydisperse = .false.
 
        thermal = dflt_int
        r0ref = dflt_real
        nb = dflt_int
 
        eu = dflt_real
        ca = dflt_real
        re_inv = dflt_real
        web = dflt_real
        poly_sigma = dflt_real
        surface_tension = .false.
 
        adv_n = .false.
 
        qbmm = .false.
        nmom = 1
        sigr = dflt_real
        sigv = dflt_real
        rhorv = 0._wp
        dist_type = dflt_int
 
        r_g = dflt_real
        r_v = dflt_real
        phi_vg = dflt_real
        phi_gv = dflt_real
        pe_c = dflt_real
        tw = dflt_real
 
        ! surface tension modeling
        sigma = dflt_real
        pi_fac = 1._wp
 
        ! Immersed Boundaries
        ib = .false.
        num_ibs = dflt_int
 
        do i = 1, num_ib_patches_max
            patch_ib(i)%geometry = dflt_int
            patch_ib(i)%x_centroid = dflt_real
            patch_ib(i)%y_centroid = dflt_real
            patch_ib(i)%z_centroid = dflt_real
            patch_ib(i)%length_x = dflt_real
            patch_ib(i)%length_y = dflt_real
            patch_ib(i)%length_z = dflt_real
            patch_ib(i)%radius = dflt_real
            patch_ib(i)%theta = dflt_real
            patch_ib(i)%c = dflt_real
            patch_ib(i)%t = dflt_real
            patch_ib(i)%m = dflt_real
            patch_ib(i)%p = dflt_real
            patch_ib(i)%slip = .false.
 
            ! Proper default values for translating STL models
            patch_ib(i)%model_scale(:) = 1._wp
            patch_ib(i)%model_translate(:) = 0._wp
            patch_ib(i)%model_rotate(:) = 0._wp
            patch_ib(i)%model_filepath(:) = dflt_char
            patch_ib(i)%model_spc = num_ray
            patch_ib(i)%model_threshold = ray_tracing_threshold
 
            ! Variables to handle moving immersed boundaries, defaulting to no movement
            patch_ib(i)%moving_ibm = 0
            patch_ib(i)%vel(:) = 0._wp
            patch_ib(i)%angles(:) = 0._wp
            patch_ib(i)%angular_vel(:) = 0._wp
            patch_ib(i)%mass = dflt_real
            patch_ib(i)%moment = dflt_real
            patch_ib(i)%centroid_offset(:) = 0._wp
 
            ! sets values of a rotation matrix which can be used when calculating rotations
            patch_ib(i)%rotation_matrix = 0._wp
            patch_ib(i)%rotation_matrix(1, 1) = 1._wp
            patch_ib(i)%rotation_matrix(2, 2) = 1._wp
            patch_ib(i)%rotation_matrix(3, 3) = 1._wp
            patch_ib(i)%rotation_matrix_inverse = patch_ib(i)%rotation_matrix
        end do
 
        chem_params%gamma_method = 1
        chem_params%transport_model = 1
 
        ! Fluids physical parameters
        do i = 1, num_fluids_max
            fluid_pp(i)%gamma = dflt_real
            fluid_pp(i)%pi_inf = dflt_real
            fluid_pp(i)%cv = 0._wp
            fluid_pp(i)%qv = 0._wp
            fluid_pp(i)%qvp = 0._wp
            fluid_pp(i)%G = 0._wp
        end do
 
        bx0 = dflt_real
 
        ! Subgrid bubble parameters
        bub_pp%R0ref = dflt_real; r0ref = dflt_real
        bub_pp%p0ref = dflt_real; p0ref = dflt_real
        bub_pp%rho0ref = dflt_real; rho0ref = dflt_real
        bub_pp%T0ref = dflt_real; t0ref = dflt_real
        bub_pp%ss = dflt_real; ss = dflt_real
        bub_pp%pv = dflt_real; pv = dflt_real
        bub_pp%vd = dflt_real; vd = dflt_real
        bub_pp%mu_l = dflt_real; mu_l = dflt_real
        bub_pp%mu_v = dflt_real; mu_v = dflt_real
        bub_pp%mu_g = dflt_real; mu_g = dflt_real
        bub_pp%gam_v = dflt_real; gam_v = dflt_real
        bub_pp%gam_g = dflt_real; gam_g = dflt_real
        bub_pp%M_v = dflt_real; m_v = dflt_real
        bub_pp%M_g = dflt_real; m_g = dflt_real
        bub_pp%k_v = dflt_real
        bub_pp%k_g = dflt_real
        bub_pp%cp_v = dflt_real; cp_v = dflt_real
        bub_pp%cp_g = dflt_real; cp_g = dflt_real
        bub_pp%R_v = dflt_real; r_v = dflt_real
        bub_pp%R_g = dflt_real; r_g = dflt_real
 
    end subroutine s_assign_default_values_to_user_inputs
 
    !> Computation of parameters, allocation procedures, and/or any other tasks needed to properly setup the module
    impure subroutine s_initialize_global_parameters_module
 
        integer :: i, j, fac
 
        if (recon_type == weno_type) then
            weno_polyn = (weno_order - 1)/2
        else if (recon_type == muscl_type) then
            muscl_polyn = muscl_order
        end if
 
        ! Determining the layout of the state vectors and overall size of the system of equations, given the dimensionality and
        ! choice of the equations of motion
 
        ! Gamma/Pi_inf Model
        if (model_eqns == 1) then
            ! Setting number of fluids
            num_fluids = 1
 
            ! Annotating structure of the state and flux vectors belonging to the system of equations defined by the selected number
            ! of spatial dimensions and the gamma/pi_inf model
            eqn_idx%cont%beg = 1
            eqn_idx%cont%end = eqn_idx%cont%beg
            eqn_idx%mom%beg = eqn_idx%cont%end + 1
            eqn_idx%mom%end = eqn_idx%cont%end + num_vels
            eqn_idx%E = eqn_idx%mom%end + 1
            eqn_idx%adv%beg = eqn_idx%E + 1
            eqn_idx%adv%end = eqn_idx%adv%beg + 1
            eqn_idx%gamma = eqn_idx%adv%beg
            eqn_idx%pi_inf = eqn_idx%adv%end
            sys_size = eqn_idx%adv%end
 
            ! Volume Fraction Model (5-equation model)
        else if (model_eqns == 2) then
            ! Annotating structure of the state and flux vectors belonging to the system of equations defined by the selected number
            ! of spatial dimensions and the volume fraction model
            eqn_idx%cont%beg = 1
            eqn_idx%cont%end = num_fluids
            eqn_idx%mom%beg = eqn_idx%cont%end + 1
            eqn_idx%mom%end = eqn_idx%cont%end + num_vels
            eqn_idx%E = eqn_idx%mom%end + 1
 
            if (igr) then
                ! Volume fractions are stored in the indices immediately following the energy equation. IGR tracks a total of (N-1)
                ! volume fractions for N fluids, hence the "-1" in eqn_idx%adv%end. If num_fluids = 1 then eqn_idx%adv%end <
                ! eqn_idx%adv%beg, which skips all loops over the volume fractions since there is no volume fraction to track
                eqn_idx%adv%beg = eqn_idx%E + 1
                eqn_idx%adv%end = eqn_idx%E + num_fluids - 1
            else
                ! Volume fractions are stored in the indices immediately following the energy equation. WENO/MUSCL + Riemann tracks
                ! a total of (N) volume fractions for N fluids, hence the lack of "-1" in eqn_idx%adv%end
                eqn_idx%adv%beg = eqn_idx%E + 1
                eqn_idx%adv%end = eqn_idx%E + num_fluids
            end if
 
            sys_size = eqn_idx%adv%end
 
            if (bubbles_euler) then
                eqn_idx%alf = eqn_idx%adv%end
            else
                eqn_idx%alf = 1
            end if
 
            if (bubbles_euler) then
                eqn_idx%bub%beg = sys_size + 1
                if (qbmm) then
                    if (nnode == 4) then
                        nmom = 6  !< Already set as a parameter
                    end if
                    eqn_idx%bub%end = eqn_idx%adv%end + nb*nmom
                else
                    if (.not. polytropic) then
                        eqn_idx%bub%end = sys_size + 4*nb
                    else
                        eqn_idx%bub%end = sys_size + 2*nb
                    end if
                end if
                sys_size = eqn_idx%bub%end
 
                if (adv_n) then
                    eqn_idx%n = eqn_idx%bub%end + 1
                    sys_size = eqn_idx%n
                end if
 
                allocate (qbmm_idx%rs(nb), qbmm_idx%vs(nb))
                allocate (qbmm_idx%ps(nb), qbmm_idx%ms(nb))
 
                if (qbmm) then
                    allocate (qbmm_idx%moms(nb, nmom))
                    allocate (qbmm_idx%fullmom(nb,0:nmom,0:nmom))
 
                    do i = 1, nb
                        do j = 1, nmom
                            qbmm_idx%moms(i, j) = eqn_idx%bub%beg + (j - 1) + (i - 1)*nmom
                        end do
                        qbmm_idx%fullmom(i, 0, 0) = qbmm_idx%moms(i, 1)
                        qbmm_idx%fullmom(i, 1, 0) = qbmm_idx%moms(i, 2)
                        qbmm_idx%fullmom(i, 0, 1) = qbmm_idx%moms(i, 3)
                        qbmm_idx%fullmom(i, 2, 0) = qbmm_idx%moms(i, 4)
                        qbmm_idx%fullmom(i, 1, 1) = qbmm_idx%moms(i, 5)
                        qbmm_idx%fullmom(i, 0, 2) = qbmm_idx%moms(i, 6)
                        qbmm_idx%rs(i) = qbmm_idx%fullmom(i, 1, 0)
                    end do
                else
                    do i = 1, nb
                        if (.not. polytropic) then
                            fac = 4
                        else
                            fac = 2
                        end if
 
                        qbmm_idx%rs(i) = eqn_idx%bub%beg + (i - 1)*fac
                        qbmm_idx%vs(i) = qbmm_idx%rs(i) + 1
 
                        if (.not. polytropic) then
                            qbmm_idx%ps(i) = qbmm_idx%vs(i) + 1
                            qbmm_idx%ms(i) = qbmm_idx%ps(i) + 1
                        end if
                    end do
                end if
            end if
 
            if (mhd) then
                eqn_idx%B%beg = sys_size + 1
                if (n == 0) then
                    eqn_idx%B%end = sys_size + 2  ! 1D: By, Bz
                else
                    eqn_idx%B%end = sys_size + 3  ! 2D/3D: Bx, By, Bz
                end if
                sys_size = eqn_idx%B%end
            end if
 
            ! Volume Fraction Model (6-equation model)
        else if (model_eqns == 3) then
            ! Annotating structure of the state and flux vectors belonging to the system of equations defined by the selected number
            ! of spatial dimensions and the volume fraction model
            eqn_idx%cont%beg = 1
            eqn_idx%cont%end = num_fluids
            eqn_idx%mom%beg = eqn_idx%cont%end + 1
            eqn_idx%mom%end = eqn_idx%cont%end + num_vels
            eqn_idx%E = eqn_idx%mom%end + 1
            eqn_idx%adv%beg = eqn_idx%E + 1
            eqn_idx%adv%end = eqn_idx%E + num_fluids
            eqn_idx%int_en%beg = eqn_idx%adv%end + 1
            eqn_idx%int_en%end = eqn_idx%adv%end + num_fluids
            sys_size = eqn_idx%int_en%end
        else if (model_eqns == 4) then
            ! 4 equation model with subgrid bubbles_euler
            eqn_idx%cont%beg = 1  ! one continuity equation
            eqn_idx%cont%end = 1  ! num_fluids
            eqn_idx%mom%beg = eqn_idx%cont%end + 1  ! one momentum equation in each direction
            eqn_idx%mom%end = eqn_idx%cont%end + num_vels
            eqn_idx%E = eqn_idx%mom%end + 1  ! one energy equation
            eqn_idx%adv%beg = eqn_idx%E + 1
            eqn_idx%adv%end = eqn_idx%adv%beg  ! one volume advection equation
            eqn_idx%alf = eqn_idx%adv%end
            sys_size = eqn_idx%alf  ! eqn_idx%adv%end
 
            if (bubbles_euler) then
                eqn_idx%bub%beg = sys_size + 1
                eqn_idx%bub%end = sys_size + 2*nb
                if (.not. polytropic) then
                    eqn_idx%bub%end = sys_size + 4*nb
                end if
                sys_size = eqn_idx%bub%end
 
                allocate (qbmm_idx%rs(nb), qbmm_idx%vs(nb))
                allocate (qbmm_idx%ps(nb), qbmm_idx%ms(nb))
                allocate (weight(nb), r0(nb))
 
                do i = 1, nb
                    if (.not. polytropic) then
                        fac = 4
                    else
                        fac = 2
                    end if
 
                    qbmm_idx%rs(i) = eqn_idx%bub%beg + (i - 1)*fac
                    qbmm_idx%vs(i) = qbmm_idx%rs(i) + 1
 
                    if (.not. polytropic) then
                        qbmm_idx%ps(i) = qbmm_idx%vs(i) + 1
                        qbmm_idx%ms(i) = qbmm_idx%ps(i) + 1
                    end if
                end do
 
                if (nb == 1) then
                    weight(:) = 1._wp
                    r0(:) = 1._wp
                else if (nb < 1) then
                    stop 'Invalid value of nb'
                end if
 
                if (polytropic) then
                    rhoref = 1._wp
                    pref = 1._wp
                end if
            end if
        end if
 
        if (model_eqns == 2 .or. model_eqns == 3) then
            if (hypoelasticity .or. hyperelasticity) then
                elasticity = .true.
                eqn_idx%stress%beg = sys_size + 1
                eqn_idx%stress%end = sys_size + (num_dims*(num_dims + 1))/2
                if (cyl_coord) eqn_idx%stress%end = eqn_idx%stress%end + 1
                ! number of stresses is 1 in 1D, 3 in 2D, 4 in 2D-Axisym, 6 in 3D
                sys_size = eqn_idx%stress%end
 
                ! shear stress index is 2 for 2D and 2,4,5 for 3D
                if (num_dims == 1) then
                    shear_num = 0
                else if (num_dims == 2) then
                    shear_num = 1
                    shear_indices(1) = eqn_idx%stress%beg - 1 + 2
                    shear_bc_flip_num = 1
                    shear_bc_flip_indices(1:2,1) = shear_indices(1)
                    ! Both x-dir and y-dir: flip tau_xy only
                else if (num_dims == 3) then
                    shear_num = 3
                    shear_indices(1:3) = eqn_idx%stress%beg - 1 + (/2, 4, 5/)
                    shear_bc_flip_num = 2
                    shear_bc_flip_indices(1,1:2) = shear_indices((/1, 2/))
                    shear_bc_flip_indices(2,1:2) = shear_indices((/1, 3/))
                    shear_bc_flip_indices(3,1:2) = shear_indices((/2, 3/))
                    ! x-dir: flip tau_xy and tau_xz y-dir: flip tau_xy and tau_yz z-dir: flip tau_xz and tau_yz
                end if
            end if
 
            if (hyperelasticity) then
                ! number of entries in the symmetric btensor plus the jacobian
                b_size = (num_dims*(num_dims + 1))/2 + 1
                tensor_size = num_dims**2 + 1
                eqn_idx%xi%beg = sys_size + 1
                eqn_idx%xi%end = sys_size + num_dims
                ! adding three more equations for the \xi field and the elastic energy
                sys_size = eqn_idx%xi%end + 1
            end if
 
            if (surface_tension) then
                eqn_idx%c = sys_size + 1
                sys_size = eqn_idx%c
            end if
 
            if (cont_damage) then
                eqn_idx%damage = sys_size + 1
                sys_size = eqn_idx%damage
            end if
 
            if (hyper_cleaning) then
                eqn_idx%psi = sys_size + 1
                sys_size = eqn_idx%psi
            end if
        end if
 
        if (chemistry) then
            eqn_idx%species%beg = sys_size + 1
            eqn_idx%species%end = sys_size + num_species
            sys_size = eqn_idx%species%end
        end if
 
        call s_configure_coordinate_bounds(recon_type, weno_polyn, muscl_polyn, igr_order, buff_size, idwint, idwbuff, viscous, &
                                           & bubbles_lagrange, m, n, p, num_dims, igr, ib)
 
#ifdef MFC_MPI
        if (qbmm .and. .not. polytropic) then
            allocate (mpi_io_data%view(1:sys_size + 2*nb*nnode))
            allocate (mpi_io_data%var(1:sys_size + 2*nb*nnode))
        else
            allocate (mpi_io_data%view(1:sys_size))
            allocate (mpi_io_data%var(1:sys_size))
        end if
 
        if (.not. down_sample) then
            do i = 1, sys_size
                allocate (mpi_io_data%var(i)%sf(0:m,0:n,0:p))
                mpi_io_data%var(i)%sf => null()
            end do
        end if
        if (qbmm .and. .not. polytropic) then
            do i = sys_size + 1, sys_size + 2*nb*nnode
                allocate (mpi_io_data%var(i)%sf(0:m,0:n,0:p))
                mpi_io_data%var(i)%sf => null()
            end do
        end if
#endif
 
        ! Allocating grid variables for the x-direction
        allocate (x_cc(0:m), x_cb(-1:m))
        ! Allocating grid variables for the y- and z-directions
        if (n > 0) then
            allocate (y_cc(0:n), y_cb(-1:n))
            if (p > 0) then
                allocate (z_cc(0:p), z_cb(-1:p))
            end if
        end if
 
        if (cyl_coord .neqv. .true.) then  ! Cartesian grid
            grid_geometry = 1
        else if (cyl_coord .and. p == 0) then  ! Axisymmetric cylindrical grid
            grid_geometry = 2
        else  ! Fully 3D cylindrical grid
            grid_geometry = 3
        end if
 
        if (.not. igr) then
            allocate (logic_grid(0:m,0:n,0:p))
        end if
 
    end subroutine s_initialize_global_parameters_module
 
    !> Configure MPI parallel I/O settings and allocate processor coordinate arrays.
    impure subroutine s_initialize_parallel_io
 
#ifdef MFC_MPI
        integer :: ierr  !< Generic flag used to identify and report MPI errors
#endif
 
        num_dims = 1 + min(1, n) + min(1, p)
 
        if (mhd) then
            num_vels = 3
        else
            num_vels = num_dims
        end if
 
        allocate (proc_coords(1:num_dims))
 
        if (parallel_io .neqv. .true.) return
 
#ifdef MFC_MPI
        ! Option for Lustre file system (Darter/Comet/Stampede)
        write (mpiiofs, '(A)') '/lustre_'
        mpiiofs = trim(mpiiofs)
        call mpi_info_create(mpi_info_int, ierr)
        call mpi_info_set(mpi_info_int, 'romio_ds_write', 'disable', ierr)
 
        ! Option for UNIX file system (Hooke/Thomson) WRITE(mpiiofs, '(A)') '/ufs_' mpiiofs = TRIM(mpiiofs) mpi_info_int =
        ! MPI_INFO_NULL
 
        allocate (start_idx(1:num_dims))
#endif
 
    end subroutine s_initialize_parallel_io
 
    !> Deallocate all global grid, index, and equation-of-state parameter arrays.
    impure subroutine s_finalize_global_parameters_module
 
        integer :: i
 
        if (bubbles_euler) then
            deallocate (qbmm_idx%rs, qbmm_idx%vs, qbmm_idx%ps, qbmm_idx%ms)
            if (qbmm) deallocate (qbmm_idx%moms, qbmm_idx%fullmom)
        end if
 
        ! Deallocating grid variables for the x-direction
 
        deallocate (x_cc, x_cb)
        ! Deallocating grid variables for the y- and z-directions
        if (n > 0) then
            deallocate (y_cc, y_cb)
            if (p > 0) then
                deallocate (z_cc, z_cb)
            end if
        end if
 
        deallocate (proc_coords)
 
#ifdef MFC_MPI
        if (parallel_io) then
            deallocate (start_idx)
            do i = 1, sys_size
                mpi_io_data%var(i)%sf => null()
            end do
 
            deallocate (mpi_io_data%var)
            deallocate (mpi_io_data%view)
        end if
#endif
 
    end subroutine s_finalize_global_parameters_module
 
end module m_global_parameters