diff --git a/src/common/m_boundary_common.fpp b/src/common/m_boundary_common.fpp index 3dc4da2c2e..10495071f7 100644 --- a/src/common/m_boundary_common.fpp +++ b/src/common/m_boundary_common.fpp @@ -1021,9 +1021,9 @@ contains integer, intent(in) :: bc_dir, bc_loc integer, intent(in) :: k, l +#ifdef MFC_SIMULATION integer :: j, i -#ifdef MFC_SIMULATION if (bc_dir == 1) then !< x-direction if (bc_loc == -1) then !bc_x%beg do i = 1, sys_size @@ -1665,7 +1665,7 @@ contains character(LEN=*), intent(in) :: step_dirpath - integer :: dir, loc, i + integer :: dir, loc character(len=path_len) :: file_path character(len=10) :: status @@ -1706,12 +1706,10 @@ contains integer :: dir, loc character(len=path_len) :: file_loc, file_path - character(len=10) :: status - #ifdef MFC_MPI integer :: ierr integer :: file_id - integer :: offset + integer(KIND=MPI_ADDRESS_KIND) :: offset character(len=7) :: proc_rank_str logical :: dir_check @@ -1770,8 +1768,6 @@ contains logical :: file_exist character(len=path_len) :: file_path - character(len=10) :: status - ! Read bc_types file_path = trim(step_dirpath)//'/bc_type.dat' inquire (FILE=trim(file_path), EXIST=file_exist) @@ -1813,12 +1809,10 @@ contains integer :: dir, loc character(len=path_len) :: file_loc, file_path - character(len=10) :: status - #ifdef MFC_MPI integer :: ierr integer :: file_id - integer :: offset + integer(KIND=MPI_ADDRESS_KIND) :: offset character(len=7) :: proc_rank_str logical :: dir_check @@ -1841,7 +1835,7 @@ contains file_path = trim(file_loc)//'/bc_'//trim(proc_rank_str)//'.dat' call MPI_File_open(MPI_COMM_SELF, trim(file_path), MPI_MODE_RDONLY, MPI_INFO_NULL, file_id, ierr) - offset = 0 + offset = int(0, KIND=MPI_ADDRESS_KIND) ! Read bc_types do dir = 1, num_dims @@ -1933,9 +1927,9 @@ contains !! boundary locations and cell-width distributions, based on !! the boundary conditions. subroutine s_populate_grid_variables_buffers - +#ifndef MFC_PRE_PROCESS integer :: i !< Generic loop iterator - +#endif #ifdef MFC_SIMULATION ! Required for compatibility between codes type(int_bounds_info) :: offset_x, offset_y, offset_z diff --git a/src/common/m_checker_common.fpp b/src/common/m_checker_common.fpp index 1ec314536a..71f0d2f442 100644 --- a/src/common/m_checker_common.fpp +++ b/src/common/m_checker_common.fpp @@ -344,8 +344,9 @@ contains !> Checks constraints on the surface tension parameters. !! Called by s_check_inputs_common for all three stages impure subroutine s_check_inputs_surface_tension - +#ifdef MFC_PRE_PROCESS integer :: i +#endif @:PROHIBIT(surface_tension .and. sigma < 0._wp, & "sigma must be greater than or equal to zero") @@ -367,7 +368,7 @@ contains @:PROHIBIT(surface_tension .and. f_is_default(patch_icpp(i)%cf_val), & "patch_icpp(i)%cf_val must be set if surface_tension is enabled") end do -#endif MFC_PRE_PROCESS +#endif end subroutine s_check_inputs_surface_tension diff --git a/src/common/m_finite_differences.fpp b/src/common/m_finite_differences.fpp index 2430374f4f..c09a47a1f8 100644 --- a/src/common/m_finite_differences.fpp +++ b/src/common/m_finite_differences.fpp @@ -70,18 +70,20 @@ contains !! @param s_cc Locations of the cell-centers in the s-coordinate direction !! @param fd_coeff_s Finite-diff. coefficients in the s-coordinate direction pure subroutine s_compute_finite_difference_coefficients(q, s_cc, fd_coeff_s, local_buff_size, & - fd_number_in, fd_order_in, offset_s) + fd_order_in, fd_number_in, offset_s) integer :: lB, lE !< loop bounds integer, intent(IN) :: q - integer, intent(IN) :: local_buff_size, fd_number_in, fd_order_in + integer, intent(IN) :: local_buff_size, fd_order_in + integer, optional, intent(IN) :: fd_number_in + type(int_bounds_info), optional, intent(IN) :: offset_s real(wp), allocatable, dimension(:, :), intent(INOUT) :: fd_coeff_s real(wp), & dimension(-local_buff_size:q + local_buff_size), & intent(IN) :: s_cc - + integer :: fd_number integer :: i !< Generic loop iterator if (present(offset_s)) then @@ -91,10 +93,15 @@ contains lB = 0 lE = q end if + if (present(fd_number_in)) then + fd_number = fd_number_in + else + fd_number = 2 + end if #ifdef MFC_POST_PROCESS if (allocated(fd_coeff_s)) deallocate (fd_coeff_s) - allocate (fd_coeff_s(-fd_number_in:fd_number_in, lb:lE)) + allocate (fd_coeff_s(-fd_number:fd_number, lb:lE)) #endif ! Computing the 1st order finite-difference coefficients diff --git a/src/common/m_helper_basic.fpp b/src/common/m_helper_basic.fpp index a1e933861a..86c30dce3a 100644 --- a/src/common/m_helper_basic.fpp +++ b/src/common/m_helper_basic.fpp @@ -7,6 +7,7 @@ module m_helper_basic use m_derived_types !< Definitions of the derived types + use m_precision_select !< Definitions of the precision types implicit none @@ -24,7 +25,7 @@ contains !> This procedure checks if two floating point numbers of wp are within tolerance. !! @param a First number. !! @param b Second number. - !! @param tol_input Relative error (default = 1.e-10_wp). + !! @param tol_input Relative error (default = 1.e-10_wp for double and 1e-6 for single). !! @return Result of the comparison. logical pure elemental function f_approx_equal(a, b, tol_input) result(res) $:GPU_ROUTINE(parallelism='[seq]') @@ -35,7 +36,11 @@ contains if (present(tol_input)) then tol = tol_input else - tol = 1.e-10_wp + if (wp == single_precision) then + tol = 1.e-6_wp + else + tol = 1.e-10_wp + end if end if if (a == b) then @@ -50,7 +55,7 @@ contains !> This procedure checks if the point numbers of wp belongs to another array are within tolerance. !! @param a First number. !! @param b Array that contains several point numbers. - !! @param tol_input Relative error (default = 1e-10_wp). + !! @param tol_input Relative error (default = 1.e-10_wp for double and 1e-6 for single). !! @return Result of the comparison. logical pure function f_approx_in_array(a, b, tol_input) result(res) $:GPU_ROUTINE(parallelism='[seq]') @@ -65,7 +70,11 @@ contains if (present(tol_input)) then tol = tol_input else - tol = 1e-10_wp + if (wp == single_precision) then + tol = 1.e-6_wp + else + tol = 1.e-10_wp + end if end if do i = 1, size(b) diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index 100c055d8d..36979deaa5 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -134,12 +134,17 @@ contains type(scalar_field), intent(in), optional :: beta integer, dimension(num_dims) :: sizes_glb, sizes_loc - integer, dimension(1) :: airfoil_glb, airfoil_loc, airfoil_start +#ifndef MFC_POST_PROCESS + integer, dimension(1) :: airfoil_glb, airfoil_loc, airfoil_start +#endif #ifdef MFC_MPI ! Generic loop iterator - integer :: i, j + integer :: i +#ifndef MFC_POST_PROCESS + integer :: j +#endif integer :: ierr !< Generic flag used to identify and report MPI errors !Altered system size for the lagrangian subgrid bubble model @@ -369,11 +374,18 @@ contains real(wp), intent(out) :: icfl_max_glb real(wp), intent(out) :: vcfl_max_glb real(wp), intent(out) :: Rc_min_glb - #ifdef MFC_SIMULATION #ifdef MFC_MPI integer :: ierr !< Generic flag used to identify and report MPI errors +#endif +#endif + ! Initiate the global variables to the local values + icfl_max_glb = icfl_max_loc + vcfl_max_glb = vcfl_max_loc + Rc_min_glb = Rc_min_loc +#ifdef MFC_SIMULATION +#ifdef MFC_MPI ! Reducing local extrema of ICFL, VCFL, CCFL and Rc numbers to their ! global extrema and bookkeeping the results on the rank 0 processor call MPI_REDUCE(icfl_max_loc, icfl_max_glb, 1, & diff --git a/src/common/m_nvtx.f90 b/src/common/m_nvtx.f90 index ce3273751a..7ebc5652a6 100644 --- a/src/common/m_nvtx.f90 +++ b/src/common/m_nvtx.f90 @@ -55,19 +55,24 @@ end subroutine nvtxRangePop subroutine nvtxStartRange(name, id) character(kind=c_char, len=*), intent(IN) :: name - integer, intent(IN), optional :: id + integer, intent(in), optional :: id + integer :: id_color +#if defined(MFC_OpenACC) && defined(__PGI) type(nvtxEventAttributes) :: event +#endif + if (present(id)) then + id_color = col(mod(id, 7) + 1) + end if + tempName = trim(name)//c_null_char #if defined(MFC_OpenACC) && defined(__PGI) - tempName = trim(name)//c_null_char - - if (.not. present(id)) then - call nvtxRangePush(tempName) - else - event%color = col(mod(id, 7) + 1) + if (present(id)) then + event%color = id_color event%message = c_loc(tempName) call nvtxRangePushEx(event) + else + call nvtxRangePushA(tempName) end if #endif diff --git a/src/common/m_phase_change.fpp b/src/common/m_phase_change.fpp index e04242a787..3b0e770477 100644 --- a/src/common/m_phase_change.fpp +++ b/src/common/m_phase_change.fpp @@ -49,8 +49,8 @@ contains !> This subroutine should dispatch to the correct relaxation solver based !! some parameter. It replaces the procedure pointer, which CCE !! is breaking on. - impure subroutine s_relaxation_solver(q_cons_vf) - type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf + impure subroutine s_relaxation_solver() + !type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf ! This is empty because in current master the procedure pointer ! was never assigned @:ASSERT(.false., "s_relaxation_solver called but it currently does nothing") diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index b9ef279f9e..404efb56c2 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -55,9 +55,11 @@ module m_variables_conversion real(wp), allocatable, dimension(:) :: Gs integer, allocatable, dimension(:) :: bubrs + +#ifdef MFC_SIMULATION real(wp), allocatable, dimension(:, :) :: Res $:GPU_DECLARE(create='[bubrs,Gs,Res]') - +#endif integer :: is1b, is2b, is3b, is1e, is2e, is3e $:GPU_DECLARE(create='[is1b,is2b,is3b,is1e,is2e,is3e]') @@ -132,9 +134,17 @@ contains real(wp) :: e_Per_Kg, Pdyn_Per_Kg real(wp) :: T_guess real(wp), dimension(1:num_species) :: Y_rs + #:if not chemistry + integer :: s !< Generic loop iterator + #:endif - integer :: s !< Generic loop iterator + ! Initiate the variables + Y_rs(:) = rhoYks(:)/rho + e_Per_Kg = energy/rho + Pdyn_Per_Kg = dyn_p/rho + E_e = 0._wp + T_guess = T #:if not chemistry ! Depending on model_eqns and bubbles_euler, the appropriate procedure ! for computing pressure is targeted by the procedure pointer @@ -154,7 +164,6 @@ contains if (hypoelasticity .and. present(G)) then ! calculate elastic contribution to Energy - E_e = 0._wp do s = stress_idx%beg, stress_idx%end if (G > 0) then E_e = E_e + ((stress/rho)**2._wp)/(4._wp*G) @@ -175,12 +184,6 @@ contains #:else - Y_rs(:) = rhoYks(:)/rho - e_Per_Kg = energy/rho - Pdyn_Per_Kg = dyn_p/rho - - T_guess = T - call get_temperature(e_Per_Kg - Pdyn_Per_Kg, T_guess, Y_rs, .true., T) call get_pressure(rho, T, Y_rs, pres) @@ -257,7 +260,10 @@ contains real(wp), optional, dimension(2), intent(out) :: Re_K - integer :: i, q + integer :: i +#ifdef MFC_SIMULATION + integer :: q +#endif real(wp), dimension(num_fluids) :: alpha_rho_K, alpha_K ! Constraining the partial densities and the volume fractions within @@ -321,14 +327,16 @@ contains qv = fluid_pp(1)%qv end if end if - + if (present(Re_K)) then + Re_K(:) = dflt_real + end if #ifdef MFC_SIMULATION ! Computing the shear and bulk Reynolds numbers from species analogs if (viscous) then if (num_fluids == 1) then ! need to consider case with num_fluids >= 2 do i = 1, 2 - Re_K(i) = dflt_real; if (Re_size(i) > 0) Re_K(i) = 0._wp + if (Re_size(i) > 0) Re_K(i) = 0._wp do q = 1, Re_size(i) Re_K(i) = (1 - alpha_K(Re_idx(i, q)))/fluid_pp(Re_idx(i, q))%Re(i) & @@ -384,8 +392,10 @@ contains real(wp), dimension(num_fluids) :: alpha_rho_K, alpha_K !< - integer :: i, j !< Generic loop iterator - + integer :: i !< Generic loop iterator +#ifdef MFC_SIMULATION + integer :: j !< Generic loop iterator +#endif ! Computing the density, the specific heat ratio function and the ! liquid stiffness function, respectively @@ -430,11 +440,14 @@ contains pi_inf = pi_inf + alpha_K(i)*pi_infs(i) qv = qv + alpha_rho_K(i)*qvs(i) end do + if (present(Re_K)) then + Re_K(:) = dflt_real + end if #ifdef MFC_SIMULATION ! Computing the shear and bulk Reynolds numbers from species analogs do i = 1, 2 - Re_K(i) = dflt_real; if (Re_size(i) > 0) Re_K(i) = 0._wp + if (Re_size(i) > 0) Re_K(i) = 0._wp do j = 1, Re_size(i) Re_K(i) = alpha_K(Re_idx(i, j))/fluid_pp(Re_idx(i, j))%Re(i) & @@ -479,22 +492,27 @@ contains real(wp), optional, intent(out) :: G_K real(wp), optional, dimension(num_fluids), intent(in) :: G - - integer :: i, j !< Generic loop iterators real(wp) :: alpha_K_sum - + integer :: i #ifdef MFC_SIMULATION - ! Constraining the partial densities and the volume fractions within - ! their physical bounds to make sure that any mixture variables that - ! are derived from them result within the limits that are set by the - ! fluids physical parameters that make up the mixture + integer :: j !< Generic loop iterators +#endif + ! Initiate the variables rho_K = 0._wp gamma_K = 0._wp pi_inf_K = 0._wp qv_K = 0._wp - alpha_K_sum = 0._wp + do i = 1, 2 + Re_K(i) = dflt_real + end do +#ifdef MFC_SIMULATION + ! Constraining the partial densities and the volume fractions within + ! their physical bounds to make sure that any mixture variables that + ! are derived from them result within the limits that are set by the + ! fluids physical parameters that make up the mixture + if (mpp_lim) then do i = 1, num_fluids alpha_rho_K(i) = max(0._wp, alpha_rho_K(i)) @@ -526,8 +544,6 @@ contains if (viscous) then do i = 1, 2 - Re_K(i) = dflt_real - if (Re_size(i) > 0) Re_K(i) = 0._wp do j = 1, Re_size(i) @@ -553,33 +569,36 @@ contains real(wp), dimension(num_fluids), intent(in) :: alpha_K, alpha_rho_K !< !! Partial densities and volume fractions - + real(wp), dimension(num_fluids) :: alpha_K_local, alpha_rho_K_local !< real(wp), dimension(2), intent(out) :: Re_K - - integer :: i, j !< Generic loop iterators - #ifdef MFC_SIMULATION + integer :: i, j !< Generic loop iterators +#endif + ! Initiate the variables rho_K = 0._wp gamma_K = 0._wp pi_inf_K = 0._wp qv_K = 0._wp - + Re_K(:) = dflt_real + alpha_K_local(:) = alpha_K(:) + alpha_rho_K_local(:) = alpha_rho_K(:) +#ifdef MFC_SIMULATION if (mpp_lim .and. (model_eqns == 2) .and. (num_fluids > 2)) then do i = 1, num_fluids - rho_K = rho_K + alpha_rho_K(i) - gamma_K = gamma_K + alpha_K(i)*gammas(i) - pi_inf_K = pi_inf_K + alpha_K(i)*pi_infs(i) - qv_K = qv_K + alpha_rho_K(i)*qvs(i) + rho_K = rho_K + alpha_rho_K_local(i) + gamma_K = gamma_K + alpha_K_local(i)*gammas(i) + pi_inf_K = pi_inf_K + alpha_K_local(i)*pi_infs(i) + qv_K = qv_K + alpha_rho_K_local(i)*qvs(i) end do else if ((model_eqns == 2) .and. (num_fluids > 2)) then do i = 1, num_fluids - 1 - rho_K = rho_K + alpha_rho_K(i) - gamma_K = gamma_K + alpha_K(i)*gammas(i) - pi_inf_K = pi_inf_K + alpha_K(i)*pi_infs(i) - qv_K = qv_K + alpha_rho_K(i)*qvs(i) + rho_K = rho_K + alpha_rho_K_local(i) + gamma_K = gamma_K + alpha_K_local(i)*gammas(i) + pi_inf_K = pi_inf_K + alpha_K_local(i)*pi_infs(i) + qv_K = qv_K + alpha_rho_K_local(i)*qvs(i) end do else - rho_K = alpha_rho_K(1) + rho_K = alpha_rho_K_local(1) gamma_K = gammas(1) pi_inf_K = pi_infs(1) qv_K = qvs(1) @@ -589,12 +608,10 @@ contains if (num_fluids == 1) then ! need to consider case with num_fluids >= 2 do i = 1, 2 - Re_K(i) = dflt_real - if (Re_size(i) > 0) Re_K(i) = 0._wp do j = 1, Re_size(i) - Re_K(i) = (1._wp - alpha_K(Re_idx(i, j)))/Res(i, j) & + Re_K(i) = (1._wp - alpha_K_local(Re_idx(i, j)))/Res(i, j) & + Re_K(i) end do @@ -604,7 +621,6 @@ contains end if end if #endif - end subroutine s_convert_species_to_mixture_variables_bubbles_acc !> The computation of parameters, the allocation of memory, @@ -612,11 +628,11 @@ contains !! other procedures that are necessary to setup the module. impure subroutine s_initialize_variables_conversion_module - integer :: i, j - - $:GPU_ENTER_DATA(copyin='[is1b,is1e,is2b,is2e,is3b,is3e]') + integer :: i #ifdef MFC_SIMULATION + integer :: j + @:ALLOCATE(gammas (1:num_fluids)) @:ALLOCATE(gs_min (1:num_fluids)) @:ALLOCATE(pi_infs(1:num_fluids)) @@ -635,6 +651,7 @@ contains @:ALLOCATE(qvps (1:num_fluids)) @:ALLOCATE(Gs (1:num_fluids)) #endif + $:GPU_ENTER_DATA(copyin='[is1b,is1e,is2b,is2e,is3b,is3e]') do i = 1, num_fluids gammas(i) = fluid_pp(i)%gamma @@ -1185,6 +1202,8 @@ contains ! Density, specific heat ratio function, liquid stiffness function ! and dynamic pressure, as defined in the incompressible flow sense, ! respectively + +#ifndef MFC_SIMULATION real(wp) :: rho real(wp) :: gamma real(wp) :: pi_inf @@ -1211,8 +1230,6 @@ contains pres_mag = 0._wp G = 0._wp - -#ifndef MFC_SIMULATION ! Converting the primitive variables to the conservative variables do l = 0, p do k = 0, n @@ -1463,6 +1480,7 @@ contains ! Partial densities, density, velocity, pressure, energy, advection ! variables, the specific heat ratio and liquid stiffness functions, ! the shear and volume Reynolds numbers and the Weber numbers +#ifdef MFC_SIMULATION real(wp), dimension(num_fluids) :: alpha_rho_K real(wp), dimension(num_fluids) :: alpha_K real(wp) :: rho_K @@ -1479,7 +1497,7 @@ contains real(wp) :: T_K, mix_mol_weight, R_gas integer :: i, j, k, l !< Generic loop iterators - +#endif is1b = is1%beg; is1e = is1%end is2b = is2%beg; is2e = is2%end is3b = is3%beg; is3e = is3%end diff --git a/src/post_process/m_data_input.f90 b/src/post_process/m_data_input.f90 index f0b4e5da6a..b039fba8ca 100644 --- a/src/post_process/m_data_input.f90 +++ b/src/post_process/m_data_input.f90 @@ -335,19 +335,14 @@ impure subroutine s_read_parallel_data_files(t_step) integer :: ifile, ierr, data_size integer, dimension(MPI_STATUS_SIZE) :: status - integer(KIND=MPI_OFFSET_KIND) :: disp integer(KIND=MPI_OFFSET_KIND) :: m_MOK, n_MOK, p_MOK - integer(KIND=MPI_OFFSET_KIND) :: WP_MOK, var_MOK, str_MOK + integer(KIND=MPI_OFFSET_KIND) :: WP_MOK, str_MOK integer(KIND=MPI_OFFSET_KIND) :: NVARS_MOK integer(KIND=MPI_OFFSET_KIND) :: MOK character(LEN=path_len + 2*name_len) :: file_loc logical :: file_exist - character(len=10) :: t_step_string - - integer :: i - allocate (x_cb_glb(-1:m_glb)) allocate (y_cb_glb(-1:n_glb)) allocate (z_cb_glb(-1:p_glb)) @@ -522,8 +517,6 @@ end subroutine s_read_parallel_conservative_data !! any other tasks needed to properly setup the module impure subroutine s_initialize_data_input_module - integer :: i !< Generic loop iterator - ! Allocating the parts of the conservative and primitive variables ! that do not require the direct knowledge of the dimensionality of ! the simulation diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp index 4c8867225e..8ca6e14956 100644 --- a/src/post_process/m_data_output.fpp +++ b/src/post_process/m_data_output.fpp @@ -724,8 +724,8 @@ contains if (precision == ${PRECISION}$) then if (p > 0) then err = DBMKOPTLIST(2, optlist) - err = DBADDIOPT(optlist, DBOPT_LO_OFFSET, lo_offset) - err = DBADDIOPT(optlist, DBOPT_HI_OFFSET, hi_offset) + err = DBADDIAOPT(optlist, DBOPT_LO_OFFSET, 3, lo_offset) + err = DBADDIAOPT(optlist, DBOPT_HI_OFFSET, 3, hi_offset) if (grid_geometry == 3) then err = DBPUTQM(dbfile, 'rectilinear_grid', 16, & 'x', 1, 'y', 1, 'z', 1, & @@ -742,8 +742,8 @@ contains err = DBFREEOPTLIST(optlist) else err = DBMKOPTLIST(2, optlist) - err = DBADDIOPT(optlist, DBOPT_LO_OFFSET, lo_offset) - err = DBADDIOPT(optlist, DBOPT_HI_OFFSET, hi_offset) + err = DBADDIAOPT(optlist, DBOPT_LO_OFFSET, 2, lo_offset) + err = DBADDIAOPT(optlist, DBOPT_HI_OFFSET, 2, hi_offset) err = DBPUTQM(dbfile, 'rectilinear_grid', 16, & 'x', 1, 'y', 1, 'z', 1, & x_cb${SFX}$, y_cb${SFX}$, DB_F77NULL, dims, 2, & diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90 index 8e7b2958f1..ed0175a968 100644 --- a/src/post_process/m_start_up.f90 +++ b/src/post_process/m_start_up.f90 @@ -245,21 +245,21 @@ impure subroutine s_save_data(t_step, varname, pres, c, H) if (omega_wrt(2) .or. omega_wrt(3) .or. qm_wrt .or. schlieren_wrt) then call s_compute_finite_difference_coefficients(m, x_cc, & fd_coeff_x, buff_size, & - fd_number, fd_order, offset_x) + fd_order, fd_number, offset_x) end if ! Computing centered finite-difference coefficients in y-direction if (omega_wrt(1) .or. omega_wrt(3) .or. qm_wrt .or. (n > 0 .and. schlieren_wrt)) then call s_compute_finite_difference_coefficients(n, y_cc, & fd_coeff_y, buff_size, & - fd_number, fd_order, offset_y) + fd_order, fd_number, offset_y) end if ! Computing centered finite-difference coefficients in z-direction if (omega_wrt(1) .or. omega_wrt(2) .or. qm_wrt .or. (p > 0 .and. schlieren_wrt)) then call s_compute_finite_difference_coefficients(p, z_cc, & fd_coeff_z, buff_size, & - fd_number, fd_order, offset_z) + fd_order, fd_number, offset_z) end if ! Adding the partial densities to the formatted database file diff --git a/src/pre_process/include/ExtrusionHardcodedIC.fpp b/src/pre_process/include/ExtrusionHardcodedIC.fpp index 264b227f21..881e63d409 100644 --- a/src/pre_process/include/ExtrusionHardcodedIC.fpp +++ b/src/pre_process/include/ExtrusionHardcodedIC.fpp @@ -37,20 +37,19 @@ #:def HardcodedDimensionsExtrusion() integer :: xRows, yRows, nRows, iix, iiy, max_files - integer :: f, iter, ios, ios2, unit, unit2, idx, idy, index_x, index_y, jump, line_count, ycount - real(wp) :: x_len, x_step, y_len, y_step + integer :: f, iter, ios, ios2, unit, unit2, idx, idy, index_x, index_y, jump, line_count + real(wp) :: x_step, y_step real(wp) :: dummy_x, dummy_y, dummy_z, x0, y0 integer :: global_offset_x, global_offset_y ! MPI subdomain offset real(wp) :: delta_x, delta_y character(len=100), dimension(sys_size) :: fileNames ! Arrays to store all data from files - character(len=200) :: errmsg + !character(len=200) :: errmsg real(wp), allocatable :: stored_values(:, :, :) real(wp), allocatable :: x_coords(:), y_coords(:) logical :: files_loaded = .false. - real(wp) :: domain_xstart, domain_xend, domain_ystart, domain_yend + real(wp) :: domain_xstart character(len=*), parameter :: init_dir = "/home/MFC/FilesDirectory" ! For example /home/MFC/examples/1D_Shock/D/ character(len=20) :: file_num_str ! For storing the file number as a string - character(len=20) :: zeros_part ! For the trailing zeros part character(len=6), parameter :: zeros_default = "000000" ! Default zeros (can be changed) #:enddef @@ -112,7 +111,7 @@ do read (unit2, *, iostat=ios2) dummy_x, dummy_y, dummy_z if (ios2 /= 0) exit - if (dummy_x == x0 .and. dummy_y /= y0) then + if (f_approx_equal(dummy_x, x0) .and. (.not. f_approx_equal(dummy_y, y0))) then yRows = yRows + 1 else exit diff --git a/src/pre_process/m_compute_levelset.fpp b/src/pre_process/m_compute_levelset.fpp index 17f66f8d68..30c2711cd2 100644 --- a/src/pre_process/m_compute_levelset.fpp +++ b/src/pre_process/m_compute_levelset.fpp @@ -157,8 +157,6 @@ contains real(wp) :: x_centroid, y_centroid, z_centroid, lz, z_max, z_min, x_act, y_act, theta real(wp), dimension(3) :: dist_vec - real(wp) :: length_z - integer :: i, j, k, l !< Loop index variables x_centroid = patch_ib(ib_patch_id)%x_centroid diff --git a/src/pre_process/m_patches.fpp b/src/pre_process/m_patches.fpp index db09af9d1f..b39bdd0776 100644 --- a/src/pre_process/m_patches.fpp +++ b/src/pre_process/m_patches.fpp @@ -49,7 +49,7 @@ module m_patches !! is to act as a pseudo volume fraction to indicate the contribution of each !! patch toward the composition of a cell's fluid state. - real(wp) :: cart_x, cart_y, cart_z + real(wp) :: cart_y, cart_z real(wp) :: sph_phi !< !! Variables to be used to hold cell locations in Cartesian coordinates if !! 3D simulation is using cylindrical coordinates @@ -130,7 +130,7 @@ contains call s_cylinder(i, ib_markers_sf, q_prim_vf, ib) call s_cylinder_levelset(levelset, levelset_norm, i) elseif (patch_ib(i)%geometry == 11) then - call s_3D_airfoil(i, ib_markers_sf, q_prim_vf, ib) + call s_3D_airfoil(i, ib_markers_sf, ib) call s_3D_airfoil_levelset(levelset, levelset_norm, i) ! STL+IBM patch elseif (patch_ib(i)%geometry == 12) then @@ -198,7 +198,7 @@ contains call s_rectangle(i, ib_markers_sf, q_prim_vf, ib) call s_rectangle_levelset(levelset, levelset_norm, i) elseif (patch_ib(i)%geometry == 4) then - call s_airfoil(i, ib_markers_sf, q_prim_vf, ib) + call s_airfoil(i, ib_markers_sf, ib) call s_airfoil_levelset(levelset, levelset_norm, i) ! STL+IBM patch elseif (patch_ib(i)%geometry == 5) then @@ -464,11 +464,10 @@ contains !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables !! @param ib True if this patch is an immersed boundary - subroutine s_airfoil(patch_id, patch_id_fp, q_prim_vf, ib_flag) + subroutine s_airfoil(patch_id, patch_id_fp, ib_flag) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp - type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf logical, optional, intent(in) :: ib_flag real(wp) :: x0, y0, f, x_act, y_act, ca_in, pa, ma, ta, theta @@ -627,11 +626,10 @@ contains !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables !! @param ib True if this patch is an immersed boundary - subroutine s_3D_airfoil(patch_id, patch_id_fp, q_prim_vf, ib_flag) + subroutine s_3D_airfoil(patch_id, patch_id_fp, ib_flag) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp - type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf logical, optional, intent(in) :: ib_flag real(wp) :: x0, y0, z0, lz, z_max, z_min, f, x_act, y_act, ca_in, pa, ma, ta, theta, xa, yt, xu, yu, xl, yl, xc, yc, dycdxc, sin_c, cos_c diff --git a/src/simulation/m_compute_cbc.fpp b/src/simulation/m_compute_cbc.fpp index 694f6735b2..71c820e871 100644 --- a/src/simulation/m_compute_cbc.fpp +++ b/src/simulation/m_compute_cbc.fpp @@ -92,7 +92,6 @@ contains real(wp), dimension(sys_size), intent(inout) :: L real(wp), intent(in) :: rho, c, dpres_ds real(wp), dimension(num_dims), intent(in) :: dvel_ds - integer :: i L(1) = f_base_L1(lambda, rho, c, dpres_ds, dvel_ds) L(2:advxe - 1) = 0._wp diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 117429487b..e087786f10 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -52,23 +52,21 @@ module m_data_output real(wp), allocatable, dimension(:, :, :) :: icfl_sf !< ICFL stability criterion real(wp), allocatable, dimension(:, :, :) :: vcfl_sf !< VCFL stability criterion - real(wp), allocatable, dimension(:, :, :) :: ccfl_sf !< CCFL stability criterion real(wp), allocatable, dimension(:, :, :) :: Rc_sf !< Rc stability criterion real(wp), public, allocatable, dimension(:, :) :: c_mass - $:GPU_DECLARE(create='[icfl_sf,vcfl_sf,ccfl_sf,Rc_sf,c_mass]') + $:GPU_DECLARE(create='[icfl_sf,vcfl_sf,Rc_sf,c_mass]') real(wp) :: icfl_max_loc, icfl_max_glb !< ICFL stability extrema on local and global grids real(wp) :: vcfl_max_loc, vcfl_max_glb !< VCFL stability extrema on local and global grids - real(wp) :: ccfl_max_loc, ccfl_max_glb !< CCFL stability extrema on local and global grids real(wp) :: Rc_min_loc, Rc_min_glb !< Rc stability extrema on local and global grids $:GPU_DECLARE(create='[icfl_max_loc,icfl_max_glb,vcfl_max_loc,vcfl_max_glb]') - $:GPU_DECLARE(create='[ccfl_max_loc,ccfl_max_glb,Rc_min_loc,Rc_min_glb]') + $:GPU_DECLARE(create='[Rc_min_loc,Rc_min_glb]') !> @name ICFL, VCFL, CCFL and Rc stability criteria extrema over all the time-steps !> @{ real(wp) :: icfl_max !< ICFL criterion maximum real(wp) :: vcfl_max !< VCFL criterion maximum - real(wp) :: ccfl_max !< CCFL criterion maximum + !real(wp) :: ccfl_max !< CCFL criterion maximum real(wp) :: Rc_min !< Rc criterion maximum !> @} @@ -215,9 +213,15 @@ contains character(LEN=path_len + 3*name_len) :: file_path !< !! Relative path to the probe data file in the case directory + character(LEN=path_len + 3*name_len) :: dir_path !< integer :: i !< Generic loop iterator logical :: file_exist + logical :: dir_exist + + write (dir_path, '(A)') trim(case_dir)//'/D' + inquire (FILE=trim(dir_path), EXIST=dir_exist) + if (.not. dir_exist) call s_create_directory(trim(dir_path)) do i = 1, num_probes ! Generating the relative path to the data file diff --git a/src/simulation/m_derived_variables.fpp b/src/simulation/m_derived_variables.fpp index 8da88a3a91..cb8897554a 100644 --- a/src/simulation/m_derived_variables.fpp +++ b/src/simulation/m_derived_variables.fpp @@ -101,17 +101,16 @@ contains end if ! Computing centered finite difference coefficients call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x, buff_size, & - fd_number, fd_order) + fd_order, fd_number) $:GPU_UPDATE(device='[fd_coeff_x]') - if (n > 0) then call s_compute_finite_difference_coefficients(n, y_cc, fd_coeff_y, buff_size, & - fd_number, fd_order) + fd_order, fd_number) $:GPU_UPDATE(device='[fd_coeff_y]') end if if (p > 0) then call s_compute_finite_difference_coefficients(p, z_cc, fd_coeff_z, buff_size, & - fd_number, fd_order) + fd_order, fd_number) $:GPU_UPDATE(device='[fd_coeff_z]') end if end if diff --git a/src/simulation/m_fftw.fpp b/src/simulation/m_fftw.fpp index 852fb90290..0d29ab0ef0 100644 --- a/src/simulation/m_fftw.fpp +++ b/src/simulation/m_fftw.fpp @@ -132,9 +132,12 @@ contains impure subroutine s_apply_fourier_filter(q_cons_vf) type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf +#if defined(MFC_OpenACC) real(c_double), pointer :: p_real(:) complex(c_double_complex), pointer :: p_cmplx(:), p_fltr_cmplx(:) - integer :: i, j, k, l !< Generic loop iterators + integer :: l !< Generic loop iterators +#endif + integer :: i, j, k !< Generic loop iterators integer :: ierr !< Generic flag used to identify and report GPU errors ! Restrict filter to processors that have cells adjacent to axis diff --git a/src/simulation/m_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp index f4a24fba7a..3a2d98519a 100644 --- a/src/simulation/m_hyperelastic.fpp +++ b/src/simulation/m_hyperelastic.fpp @@ -71,16 +71,16 @@ contains ! Computing centered finite difference coefficients call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x, buff_size, & - fd_number, fd_order) + fd_order, fd_number) $:GPU_UPDATE(device='[fd_coeff_x]') if (n > 0) then call s_compute_finite_difference_coefficients(n, y_cc, fd_coeff_y, buff_size, & - fd_number, fd_order) + fd_order, fd_number) $:GPU_UPDATE(device='[fd_coeff_y]') end if if (p > 0) then call s_compute_finite_difference_coefficients(p, z_cc, fd_coeff_z, buff_size, & - fd_number, fd_order) + fd_order, fd_number) $:GPU_UPDATE(device='[fd_coeff_z]') end if diff --git a/src/simulation/m_hypoelastic.fpp b/src/simulation/m_hypoelastic.fpp index 3f736b0b0b..52087e7b43 100644 --- a/src/simulation/m_hypoelastic.fpp +++ b/src/simulation/m_hypoelastic.fpp @@ -67,16 +67,16 @@ contains ! Computing centered finite difference coefficients call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x_h, buff_size, & - fd_number, fd_order) + fd_order, fd_number) $:GPU_UPDATE(device='[fd_coeff_x_h]') if (n > 0) then call s_compute_finite_difference_coefficients(n, y_cc, fd_coeff_y_h, buff_size, & - fd_number, fd_order) + fd_order, fd_number) $:GPU_UPDATE(device='[fd_coeff_y_h]') end if if (p > 0) then call s_compute_finite_difference_coefficients(p, z_cc, fd_coeff_z_h, buff_size, & - fd_number, fd_order) + fd_order, fd_number) $:GPU_UPDATE(device='[fd_coeff_z_h]') end if diff --git a/src/simulation/m_mhd.fpp b/src/simulation/m_mhd.fpp index 8112b3af7e..93acb73a28 100644 --- a/src/simulation/m_mhd.fpp +++ b/src/simulation/m_mhd.fpp @@ -21,10 +21,10 @@ module m_mhd s_finalize_mhd_powell_module, & s_compute_mhd_powell_rhs - real(wp), allocatable, dimension(:, :, :) :: du_dx, du_dy, du_dz - real(wp), allocatable, dimension(:, :, :) :: dv_dx, dv_dy, dv_dz + real(wp), allocatable, dimension(:, :, :) :: du_dx, du_dy + real(wp), allocatable, dimension(:, :, :) :: dv_dx, dv_dy real(wp), allocatable, dimension(:, :, :) :: dw_dx, dw_dy, dw_dz - $:GPU_DECLARE(create='[du_dx,du_dy,du_dz,dv_dx,dv_dy,dv_dz,dw_dx,dw_dy,dw_dz]') + $:GPU_DECLARE(create='[du_dx,du_dy,dv_dx,dv_dy,dw_dx,dw_dy,dw_dz]') real(wp), allocatable, dimension(:, :) :: fd_coeff_x_h real(wp), allocatable, dimension(:, :) :: fd_coeff_y_h @@ -51,12 +51,12 @@ contains end if ! Computing centered finite difference coefficients - call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x_h, buff_size, fd_number, fd_order) + call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x_h, buff_size, fd_order, fd_number) $:GPU_UPDATE(device='[fd_coeff_x_h]') - call s_compute_finite_difference_coefficients(n, y_cc, fd_coeff_y_h, buff_size, fd_number, fd_order) + call s_compute_finite_difference_coefficients(n, y_cc, fd_coeff_y_h, buff_size, fd_order, fd_number) $:GPU_UPDATE(device='[fd_coeff_y_h]') if (p > 0) then - call s_compute_finite_difference_coefficients(p, z_cc, fd_coeff_z_h, buff_size, fd_number, fd_order) + call s_compute_finite_difference_coefficients(p, z_cc, fd_coeff_z_h, buff_size, fd_order, fd_number) $:GPU_UPDATE(device='[fd_coeff_z_h]') end if diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index f5cc89b4c5..8c441cfe49 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -247,7 +247,7 @@ contains integer, intent(in) :: mpi_dir, pbc_loc - integer :: i, j, k, l, r, q !< Generic loop iterators + integer :: j, k, l, r !< Generic loop iterators integer :: buffer_counts(1:3), buffer_count @@ -255,7 +255,7 @@ contains integer :: beg_end(1:2), grid_dims(1:3) integer :: dst_proc, src_proc, recv_tag, send_tag - logical :: beg_end_geq_0, qbmm_comm + logical :: beg_end_geq_0 integer :: pack_offset, unpack_offset diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index e6a90499a5..e10219187d 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -113,21 +113,6 @@ module m_rhs type(scalar_field), allocatable, dimension(:) :: tau_Re_vf $:GPU_DECLARE(create='[tau_Re_vf]') - type(vector_field) :: gm_alpha_qp !< - !! The gradient magnitude of the volume fractions at cell-interior Gaussian - !! quadrature points. gm_alpha_qp is calculated from individual first-order - !! spatial derivatives located in dq_prim_ds_qp. - - $:GPU_DECLARE(create='[gm_alpha_qp]') - - !> @name The left and right WENO-reconstructed cell-boundary values of the cell- - !! average gradient magnitude of volume fractions, located in gm_alpha_qp. - !> @{ - type(vector_field), allocatable, dimension(:) :: gm_alphaL_n - type(vector_field), allocatable, dimension(:) :: gm_alphaR_n - $:GPU_DECLARE(create='[gm_alphaL_n,gm_alphaR_n]') - !> @} - !> @name The cell-boundary values of the fluxes (src - source, gsrc - geometrical !! source). These are computed by applying the chosen Riemann problem solver !! .on the left and right cell-boundary values of the primitive variables @@ -602,10 +587,6 @@ contains end if ! END: Allocation/Association of qK_cons_n and qK_prim_n - ! Allocation of gm_alphaK_n - @:ALLOCATE(gm_alphaL_n(1:num_dims)) - @:ALLOCATE(gm_alphaR_n(1:num_dims)) - if (alt_soundspeed) then @:ALLOCATE(blkmod1(0:m, 0:n, 0:p), blkmod2(0:m, 0:n, 0:p), alpha1(0:m, 0:n, 0:p), alpha2(0:m, 0:n, 0:p), Kterm(0:m, 0:n, 0:p)) end if diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 8196403564..23d1967b43 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -312,7 +312,6 @@ contains real(wp), dimension(6) :: tau_e_L, tau_e_R real(wp) :: G_L, G_R real(wp), dimension(2) :: Re_L, Re_R - real(wp), dimension(3) :: xi_field_L, xi_field_R real(wp) :: rho_avg real(wp) :: H_avg @@ -360,7 +359,7 @@ contains $:GPU_PARALLEL_LOOP(collapse=3, private='[alpha_rho_L, alpha_rho_R, & & vel_L, vel_R, alpha_L, alpha_R, tau_e_L, tau_e_R, & & G_L, G_R, Re_L, Re_R, rho_avg, h_avg, gamma_avg, & - & s_L, s_R, s_S, Ys_L, Ys_R, xi_field_L, xi_field_R, & + & s_L, s_R, s_S, Ys_L, Ys_R, & & Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR, & & Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2, c_fast, & & pres_mag, B, Ga, vdotB, B2, b4, cm, pcorr, & @@ -3671,7 +3670,6 @@ contains !! @param ix Index bounds in the x-dir !! @param iy Index bounds in the y-dir !! @param iz Index bounds in the z-dir - !! @param q_prim_vf Cell-averaged primitive variables subroutine s_initialize_riemann_solver( & flux_src_vf, & norm_dir) diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index e004252060..0d8878abf8 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -1236,7 +1236,7 @@ contains call nvtxEndRange call cpu_time(finish) if (cfl_dt) then - nt = mytime/t_save + nt = int(mytime/t_save) else nt = int((t_step - t_step_start)/(t_step_save)) end if @@ -1336,8 +1336,9 @@ contains end subroutine s_initialize_modules impure subroutine s_initialize_mpi_domain - integer :: ierr #ifdef MFC_OpenACC + integer :: ierr + real(wp) :: starttime, endtime integer :: num_devices, local_size, num_nodes, ppn, my_device_num integer :: dev, devNum, local_rank diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index d040650bfa..833c32b568 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -945,8 +945,6 @@ contains integer, intent(in) :: stage - type(vector_field) :: gm_alpha_qp - call s_convert_conservative_to_primitive_variables( & q_cons_ts(1)%vf, & q_T_sf, & @@ -989,7 +987,6 @@ contains real(wp) :: c !< Cell-avg. sound speed real(wp) :: H !< Cell-avg. enthalpy real(wp), dimension(2) :: Re !< Cell-avg. Reynolds numbers - type(vector_field) :: gm_alpha_qp real(wp) :: dt_local integer :: j, k, l !< Generic loop iterators diff --git a/src/simulation/m_weno.fpp b/src/simulation/m_weno.fpp index 56beaea979..e6814e0598 100644 --- a/src/simulation/m_weno.fpp +++ b/src/simulation/m_weno.fpp @@ -1221,9 +1221,6 @@ contains !! Determines the amount of freedom available from utilizing a large !! value for the local curvature. The default value for beta is 4/3. - real(wp), parameter :: alpha_mp = 2._wp - real(wp), parameter :: beta_mp = 4._wp/3._wp - $:GPU_PARALLEL_LOOP(collapse=4,private='[d]') do l = is3_weno%beg, is3_weno%end do k = is2_weno%beg, is2_weno%end @@ -1256,7 +1253,7 @@ contains vL_UL = v_rs_ws(j, k, l, i) & - (v_rs_ws(j + 1, k, l, i) & - - v_rs_ws(j, k, l, i))*alpha_mp + - v_rs_ws(j, k, l, i))*alpha vL_MD = (v_rs_ws(j, k, l, i) & + v_rs_ws(j - 1, k, l, i) & @@ -1264,7 +1261,7 @@ contains vL_LC = v_rs_ws(j, k, l, i) & - (v_rs_ws(j + 1, k, l, i) & - - v_rs_ws(j, k, l, i))*5.e-1_wp + beta_mp*d_LC + - v_rs_ws(j, k, l, i))*5.e-1_wp + beta*d_LC vL_min = max(min(v_rs_ws(j, k, l, i), & v_rs_ws(j - 1, k, l, i), & @@ -1315,7 +1312,7 @@ contains vR_UL = v_rs_ws(j, k, l, i) & + (v_rs_ws(j, k, l, i) & - - v_rs_ws(j - 1, k, l, i))*alpha_mp + - v_rs_ws(j - 1, k, l, i))*alpha vR_MD = (v_rs_ws(j, k, l, i) & + v_rs_ws(j + 1, k, l, i) & @@ -1323,7 +1320,7 @@ contains vR_LC = v_rs_ws(j, k, l, i) & + (v_rs_ws(j, k, l, i) & - - v_rs_ws(j - 1, k, l, i))*5.e-1_wp + beta_mp*d_LC + - v_rs_ws(j - 1, k, l, i))*5.e-1_wp + beta*d_LC vR_min = max(min(v_rs_ws(j, k, l, i), & v_rs_ws(j + 1, k, l, i), &