From edd8d8bd8e463515d67aed59f558ed902110bcfd Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Wed, 11 Sep 2024 09:53:41 -0400 Subject: [PATCH 01/10] add comprehensive debugging --- src/simulation/m_global_parameters.fpp | 12 ++-- src/simulation/m_mpi_proxy.fpp | 2 +- src/simulation/m_sim_helpers.f90 | 2 + src/simulation/m_start_up.fpp | 3 +- src/simulation/m_time_steppers.fpp | 79 ++++++++++++++++++++++++++ toolchain/mfc/run/case_dicts.py | 1 + 6 files changed, 93 insertions(+), 6 deletions(-) diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 8797a8f45..7efb1cdb8 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -158,6 +158,7 @@ module m_global_parameters logical :: mixture_err !< Mixture properties correction logical :: hypoelasticity !< hypoelasticity modeling logical :: cu_tensor + logical :: comp_debug !< Variable checking at every RK step logical :: bodyForces logical :: bf_x, bf_y, bf_z !< body force toggle in three directions @@ -624,6 +625,9 @@ contains ! Cuda aware MPI cu_tensor = .false. + ! Comprehensive debugging + comp_debug = .false. + bodyForces = .false. bf_x = .false.; bf_y = .false.; bf_z = .false. !< amplitude, frequency, and phase shift sinusoid in each direction @@ -1088,12 +1092,12 @@ contains @:ALLOCATE_GLOBAL(x_cc(-buff_size:m + buff_size)) @:ALLOCATE_GLOBAL(dx(-buff_size:m + buff_size)) - if (n == 0) return; + if (n == 0) return; @:ALLOCATE_GLOBAL(y_cb(-1 - buff_size:n + buff_size)) @:ALLOCATE_GLOBAL(y_cc(-buff_size:n + buff_size)) @:ALLOCATE_GLOBAL(dy(-buff_size:n + buff_size)) - if (p == 0) return; + if (p == 0) return; @:ALLOCATE_GLOBAL(z_cb(-1 - buff_size:p + buff_size)) @:ALLOCATE_GLOBAL(z_cc(-buff_size:p + buff_size)) @:ALLOCATE_GLOBAL(dz(-buff_size:p + buff_size)) @@ -1159,10 +1163,10 @@ contains ! Deallocating grid variables for the x-, y- and z-directions @:DEALLOCATE_GLOBAL(x_cb, x_cc, dx) - if (n == 0) return; + if (n == 0) return; @:DEALLOCATE_GLOBAL(y_cb, y_cc, dy) - if (p == 0) return; + if (p == 0) return; @:DEALLOCATE_GLOBAL(z_cb, z_cc, dz) end subroutine s_finalize_global_parameters_module diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index b395fc624..7970ac8d3 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -203,7 +203,7 @@ contains & 'polydisperse', 'qbmm', 'acoustic_source', 'probe_wrt', 'integral_wrt', & & 'prim_vars_wrt', 'weno_avg', 'file_per_process', 'relax', & & 'adv_n', 'adap_dt', 'ib', 'bodyForces', 'bf_x', 'bf_y', 'bf_z', & - & 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt' ] + & 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt', 'comp_debug' ] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor diff --git a/src/simulation/m_sim_helpers.f90 b/src/simulation/m_sim_helpers.f90 index eca90a32e..be95197f6 100644 --- a/src/simulation/m_sim_helpers.f90 +++ b/src/simulation/m_sim_helpers.f90 @@ -6,6 +6,8 @@ module m_sim_helpers use m_global_parameters use m_variables_conversion + + use m_mpi_proxy ! ========================================================================== implicit none diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 4feb49a42..bfd09ce9f 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -162,7 +162,8 @@ contains pi_fac, adv_n, adap_dt, bf_x, bf_y, bf_z, & k_x, k_y, k_z, w_x, w_y, w_z, p_x, p_y, p_z, & g_x, g_y, g_z, n_start, t_save, t_stop, & - cfl_adap_dt, cfl_const_dt, cfl_target + cfl_adap_dt, cfl_const_dt, cfl_target, & + comp_debug ! Checking that an input file has been provided by the user. If it ! has, then the input file is read in, otherwise, simulation exits. diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index a4ff54a58..8c9d9c41a 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -38,6 +38,10 @@ module m_time_steppers use m_nvtx use m_body_forces + + use ieee_arithmetic + + use m_sim_helpers ! ========================================================================== implicit none @@ -391,6 +395,8 @@ contains if (model_eqns == 3) call s_pressure_relaxation_procedure(q_cons_ts(1)%vf) + if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, q_prim_vf, t_step, 1) + if (adv_n) call s_comp_alpha_from_n(q_cons_ts(1)%vf) if (ib) then @@ -498,6 +504,8 @@ contains call s_pressure_relaxation_procedure(q_cons_ts(2)%vf) end if + if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, q_prim_vf, t_step, 1) + if (adv_n) call s_comp_alpha_from_n(q_cons_ts(2)%vf) if (ib) then @@ -573,6 +581,8 @@ contains call s_pressure_relaxation_procedure(q_cons_ts(1)%vf) end if + if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, q_prim_vf, t_step, 2) + if (adv_n) call s_comp_alpha_from_n(q_cons_ts(1)%vf) if (ib) then @@ -683,6 +693,8 @@ contains call s_pressure_relaxation_procedure(q_cons_ts(2)%vf) end if + if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, q_prim_vf, t_step, 1) + if (adv_n) call s_comp_alpha_from_n(q_cons_ts(2)%vf) if (ib) then @@ -758,6 +770,8 @@ contains call s_pressure_relaxation_procedure(q_cons_ts(2)%vf) end if + if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, q_prim_vf, t_step, 2) + if (adv_n) call s_comp_alpha_from_n(q_cons_ts(2)%vf) if (ib) then @@ -832,6 +846,8 @@ contains call s_pressure_relaxation_procedure(q_cons_ts(1)%vf) end if + if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, q_prim_vf, t_step, 3) + if (adv_n) call s_comp_alpha_from_n(q_cons_ts(1)%vf) if (ib) then @@ -994,6 +1010,68 @@ contains end subroutine s_apply_bodyforces + subroutine s_comprehensive_debug(q_cons_vf, q_prim_vf, t_step, stage) + + type(scalar_field), dimension(sys_size) :: q_cons_vf, q_prim_vf + integer, intent(in) :: t_step, stage + integer :: j, k, l, i + integer errors + + errors = 0 + + ! Check all variables for NaNs + do i = 1, sys_size + !$acc update host(q_cons_vf(i)%sf) + do l = 0, p + do k = 0, n + do j = 0, m + if (ieee_is_nan(q_cons_vf(i)%sf(j, k, l))) then + print *, "NaN(s) in conservative variables after RK stage", stage, "at", j, k, l, i, proc_rank, t_step, m, n, p + errors = errors + 1 + end if + end do + end do + end do + end do + + ! Check for invalid volume fractions + do i = advxb, advxe + do l = 0, p + do k = 0, n + do j = 0, m + if (q_cons_vf(i)%sf(j, k, l) < 0d0) then + print *, "Volume Fraction is Negative after RK stage", stage, "at", j, k, l, i, proc_rank, t_step, m, n, p + errors = errors + 1 + elseif (q_cons_Vf(i)%sf(j, k, l) > 1d0) then + print *, "Volume Fraction is greater than 1 after RK stage", stage, "at", j, k, l, i, proc_rank, t_step, m, n, p + errors = errors + 1 + end if + end do + end do + end do + end do + + ! Check for invalid densities + do i = contxb, contxe + do l = 0, p + do k = 0, n + do j = 0, m + if (q_cons_vf(advxb + i -1)%sf(j, k, l) < 0d0 .and. q_cons_vf(i)%sf(j, k, l) > 0d0) then + print *, "Density is negative after RK stage", stage, "at", j, k, l, i, proc_rank, t_step, m, n, p + errors = errors + 1 + end if + end do + end do + end do + end do + + if (errors /= 0) then + call s_write_data_files(q_cons_vf, q_prim_vf, t_step) + call s_mpi_abort("Errors found in conservative variables") + endif + + end subroutine s_comprehensive_debug + !> This subroutine saves the temporary q_prim_vf vector !! into the q_prim_ts vector that is then used in p_main !! @param t_step current time-step @@ -1033,6 +1111,7 @@ contains end if end subroutine s_time_step_cycling + !> Module deallocation and/or disassociation procedures subroutine s_finalize_time_steppers_module diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index 95b2a5bcc..d574d41e1 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -219,6 +219,7 @@ def analytic(self): 't_save': ParamType.REAL, 'cfl_target': ParamType.REAL, 'low_Mach': ParamType.INT, + 'comp_debug': ParamType.LOG, }) # NOTE: Not currently present From 1117c5c2d6eb6e934322c10e6cbd88e5d7d511b1 Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Tue, 17 Sep 2024 14:26:30 -0400 Subject: [PATCH 02/10] comp_debug on CPUs --- src/post_process/m_global_parameters.fpp | 2 + src/post_process/m_mpi_proxy.fpp | 2 +- src/post_process/m_start_up.f90 | 52 +++++++++++++++++++++++- src/post_process/p_main.fpp | 2 + src/simulation/m_time_steppers.fpp | 39 +++++++++++++++--- toolchain/mfc/run/case_dicts.py | 1 + 6 files changed, 90 insertions(+), 8 deletions(-) diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 8a6b2805f..84172bb55 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -104,6 +104,7 @@ module m_global_parameters logical :: mixture_err !< Mixture error limiter logical :: alt_soundspeed !< Alternate sound speed logical :: hypoelasticity !< Turn hypoelasticity on + logical :: comp_debug !< Turn on comprehensive debug !> @} !> @name Annotations of the structure, i.e. the organization, of the state vectors @@ -305,6 +306,7 @@ contains relax = .false. relax_model = dflt_int hypoelasticity = .false. + comp_debug = .false. bc_x%beg = dflt_int; bc_x%end = dflt_int bc_y%beg = dflt_int; bc_y%end = dflt_int diff --git a/src/post_process/m_mpi_proxy.fpp b/src/post_process/m_mpi_proxy.fpp index 825f61b49..91b3a91ed 100644 --- a/src/post_process/m_mpi_proxy.fpp +++ b/src/post_process/m_mpi_proxy.fpp @@ -169,7 +169,7 @@ contains & 'heat_ratio_wrt', 'pi_inf_wrt', 'pres_inf_wrt', 'cons_vars_wrt', & & 'prim_vars_wrt', 'c_wrt', 'qm_wrt','schlieren_wrt', 'bubbles', 'qbmm', & & 'polytropic', 'polydisperse', 'file_per_process', 'relax', 'cf_wrt', & - & 'adv_n', 'ib', 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt' ] + & 'adv_n', 'ib', 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt', 'comp_debug' ] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90 index 2dd94da62..10016f4cf 100644 --- a/src/post_process/m_start_up.f90 +++ b/src/post_process/m_start_up.f90 @@ -76,7 +76,7 @@ subroutine s_read_input_file polydisperse, poly_sigma, file_per_process, relax, & relax_model, cf_wrt, sigma, adv_n, ib, & cfl_adap_dt, cfl_const_dt, t_save, t_stop, n_start, & - cfl_target + cfl_target, comp_debug ! Inquiring the status of the post_process.inp file file_loc = 'post_process.inp' @@ -176,6 +176,56 @@ subroutine s_perform_time_step(t_step) call s_convert_conservative_to_primitive_variables(q_cons_vf, q_prim_vf) end subroutine s_perform_time_step + subroutine s_perform_comprehensive_debug(varname, pres, c, H) + + real(kind(0d0)) :: pres + real(kind(0d0)) :: c + real(kind(0d0)) :: H + integer :: t_fail, ios + logical :: exists + character(LEN=name_len) :: varname !< + character(LEN=name_len) :: file_name = 'comp_debug.txt' + character(LEN=path_len + name_len) :: file_path + + ! Opening the run-time information file + file_path = trim(case_dir)//'/'//trim(file_name) + + inquire(file=file_path, exist=exists) + if (exists) then + open(12, file=file_path) + ! Read the file line by line + do + read(12, '(I9)', iostat=ios) t_fail + if (ios /= 0) exit ! Exit loop on error or end-of-file + end do + if (t_fail == -1) return + else + return + end if + + if (proc_rank == 0) then + print("(A, I6)"), " Post Processing suspicious time-step: ", t_fail + end if + + ! Populating the grid and conservative variables + call s_read_data_files(t_fail) + ! Populating the buffer regions of the grid variables + if (buff_size > 0) then + call s_populate_grid_variables_buffer_regions() + end if + + ! Populating the buffer regions of the conservative variables + if (buff_size > 0) then + call s_populate_conservative_variables_buffer_regions() + end if + + ! Converting the conservative variables to the primitive ones + call s_convert_conservative_to_primitive_variables(q_cons_vf, q_prim_vf) + + call s_save_data(t_fail, varname, pres, c, H) + + end subroutine s_perform_comprehensive_debug + subroutine s_save_data(t_step, varname, pres, c, H) integer, intent(inout) :: t_step diff --git a/src/post_process/p_main.fpp b/src/post_process/p_main.fpp index e48bb8011..094eb9959 100644 --- a/src/post_process/p_main.fpp +++ b/src/post_process/p_main.fpp @@ -34,6 +34,8 @@ program p_main call s_initialize_modules() + if (comp_debug) call s_perform_comprehensive_debug(varname, pres, c, H) + if (cfl_dt) then t_step = n_start n_save = int(t_stop/t_save) + 1 diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 8c9d9c41a..f9d904737 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -1016,6 +1016,19 @@ contains integer, intent(in) :: t_step, stage integer :: j, k, l, i integer errors + logical :: exists + + character(LEN=name_len) :: file_name = 'comp_debug.txt' + character(LEN=path_len + name_len) :: file_path + character(100) :: str_format + + ! Opening the run-time information file + file_path = trim(case_dir)//'/'//trim(file_name) + + str_format = "(I9, A, I3, A, I4, I4, I4, A, I2, A, I5, A, I5, I5, I5)" + + open (12, FILE=trim(file_path), & + STATUS='replace') errors = 0 @@ -1026,7 +1039,9 @@ contains do k = 0, n do j = 0, m if (ieee_is_nan(q_cons_vf(i)%sf(j, k, l))) then - print *, "NaN(s) in conservative variables after RK stage", stage, "at", j, k, l, i, proc_rank, t_step, m, n, p + write(12, str_format) t_step, " NaN(s) in conservative variables after RK stage ", & + stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & + " (m, n, p)", m, n, p errors = errors + 1 end if end do @@ -1040,10 +1055,14 @@ contains do k = 0, n do j = 0, m if (q_cons_vf(i)%sf(j, k, l) < 0d0) then - print *, "Volume Fraction is Negative after RK stage", stage, "at", j, k, l, i, proc_rank, t_step, m, n, p + write(12, str_format) t_step, " Volume fraction < 0 after RK stage ", & + stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & + " (m, n, p)", m, n, p errors = errors + 1 - elseif (q_cons_Vf(i)%sf(j, k, l) > 1d0) then - print *, "Volume Fraction is greater than 1 after RK stage", stage, "at", j, k, l, i, proc_rank, t_step, m, n, p + elseif (q_cons_vf(i)%sf(j, k, l) > 1d0 + verysmall) then + write(12, str_format) t_step, " Volume fraction > 1 after RK stage ", & + stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & + " (m, n, p)", m, n, p errors = errors + 1 end if end do @@ -1056,8 +1075,12 @@ contains do l = 0, p do k = 0, n do j = 0, m - if (q_cons_vf(advxb + i -1)%sf(j, k, l) < 0d0 .and. q_cons_vf(i)%sf(j, k, l) > 0d0) then - print *, "Density is negative after RK stage", stage, "at", j, k, l, i, proc_rank, t_step, m, n, p + if (q_cons_vf(advxb + i -1)%sf(j, k, l) < 0d0 .and. q_cons_vf(i)%sf(j, k, l) < 0d0 .or. & + q_cons_vf(advxb + i -1)%sf(j, k, l) > 0d0 .and. q_cons_Vf(i)%sf(j, k, l) < 0d0) then + print*, q_cons_vf(advxb + i - 1)%sf(j, k, l), q_cons_vf(i)%sf(j, k, l) + write(12, str_format) t_step, " Density is negative after RK stage ", & + stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & + " (m, n, p)", m, n, p errors = errors + 1 end if end do @@ -1066,10 +1089,14 @@ contains end do if (errors /= 0) then + close(12) call s_write_data_files(q_cons_vf, q_prim_vf, t_step) call s_mpi_abort("Errors found in conservative variables") endif + write(12, "(I3)") -1 + close(12) + end subroutine s_comprehensive_debug !> This subroutine saves the temporary q_prim_vf vector diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index d574d41e1..610bbbe46 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -335,6 +335,7 @@ def analytic(self): 't_save': ParamType.REAL, 't_stop': ParamType.REAL, 'n_start': ParamType.INT, + 'comp_debug': ParamType.LOG, }) for cmp_id in range(1,3+1): From 92b9a8fb45729b9178ec6f4cfdaf8d4d4e687da5 Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Tue, 17 Sep 2024 14:27:09 -0400 Subject: [PATCH 03/10] minor bug fix --- src/simulation/m_start_up.fpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index bfd09ce9f..e06685e13 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -1105,13 +1105,13 @@ contains integer :: i, j, k, l if (cfl_dt) then - if (cfl_const_dt .and. t_step == 1) call s_compute_dt() + if (cfl_const_dt .and. t_step == 0) call s_compute_dt() if (cfl_adap_dt) call s_compute_dt() if (t_step == 0) dt_init = dt - if (dt < 1d-3*dt_init) call s_mpi_abort("Delta t has become too small") + if (dt < 1d-3*dt_init .and. cfl_adap_dt) call s_mpi_abort("Delta t has become too small") end if if (cfl_dt) then From 73f7931555d8dfd6cfe729c907b8048f8508dbce Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Thu, 19 Sep 2024 17:12:29 -0400 Subject: [PATCH 04/10] documentation --- docs/documentation/case.md | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index a4e80117f..ccb405e33 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -90,11 +90,18 @@ Definition of the parameters is described in the following subsections. | ---: | :----: | :--- | | `run_time_info` | Logical | Output run-time information | | `rdma_mpi` | Logical | (GPUs) Enable RDMA for MPI communication. | +| `comp_debug` | Logical | Comprehensive variable checking | - `run_time_info` generates a text file that includes run-time information including the CFL number(s) at each time-step. - `rdma_mpi` optimizes data transfers between GPUs using Remote Direct Memory Access (RDMA). The underlying MPI implementation and communication infrastructure must support this feature, detecting GPU pointers and performing RDMA accordingly. +- `comp_debug` enables comprehensive error checking. +At each Runge-Kutta sub-step, all conservative variables are checked for NaNs. +The volume fractions are checked to ensure they are in the range [0, 1]. +Negative densities are also checked for. +If any of these checks find problems, the file `comp_debug.txt` will be written to the case directory with information +about what problems were found and the simulation state will be saved for visualization. ### 2. Computational Domain From f4283a4af8af992fd76c9806bdf1b0ee5883a66f Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Sun, 22 Sep 2024 15:25:59 -0400 Subject: [PATCH 05/10] refactor --- src/simulation/m_data_output.fpp | 291 ------- src/simulation/m_sim_helpers.f90 | 250 ------ src/simulation/m_sim_helpers.fpp | 1272 ++++++++++++++++++++++++++++ src/simulation/m_start_up.fpp | 628 +------------- src/simulation/m_time_steppers.fpp | 91 -- 5 files changed, 1278 insertions(+), 1254 deletions(-) delete mode 100644 src/simulation/m_sim_helpers.f90 create mode 100644 src/simulation/m_sim_helpers.fpp diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index b5d33da11..0a48562c1 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -26,8 +26,6 @@ module m_data_output use m_helper - use m_sim_helpers - use m_delay_file_access use m_ibm @@ -37,14 +35,11 @@ module m_data_output private; public :: s_initialize_data_output_module, & - s_open_run_time_information_file, & s_open_probe_files, & - s_write_run_time_information, & s_write_data_files, & s_write_serial_data_files, & s_write_parallel_data_files, & s_write_probe_files, & - s_close_run_time_information_file, & s_close_probe_files, & s_finalize_data_output_module @@ -70,103 +65,11 @@ module m_data_output end subroutine s_write_abstract_data_files end interface ! ======================================================== -#ifdef CRAY_ACC_WAR - @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), icfl_sf) - @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), vcfl_sf) - @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), ccfl_sf) - @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), Rc_sf) - !$acc declare link(icfl_sf, vcfl_sf, ccfl_sf, Rc_sf) -#else - real(kind(0d0)), allocatable, dimension(:, :, :) :: icfl_sf !< ICFL stability criterion - real(kind(0d0)), allocatable, dimension(:, :, :) :: vcfl_sf !< VCFL stability criterion - real(kind(0d0)), allocatable, dimension(:, :, :) :: ccfl_sf !< CCFL stability criterion - real(kind(0d0)), allocatable, dimension(:, :, :) :: Rc_sf !< Rc stability criterion - !$acc declare create(icfl_sf, vcfl_sf, ccfl_sf, Rc_sf) -#endif - - real(kind(0d0)) :: icfl_max_loc, icfl_max_glb !< ICFL stability extrema on local and global grids - real(kind(0d0)) :: vcfl_max_loc, vcfl_max_glb !< VCFL stability extrema on local and global grids - real(kind(0d0)) :: ccfl_max_loc, ccfl_max_glb !< CCFL stability extrema on local and global grids - real(kind(0d0)) :: Rc_min_loc, Rc_min_glb !< Rc stability extrema on local and global grids - !$acc declare create(icfl_max_loc, icfl_max_glb, vcfl_max_loc, vcfl_max_glb, ccfl_max_loc, ccfl_max_glb, Rc_min_loc, Rc_min_glb) - - !> @name ICFL, VCFL, CCFL and Rc stability criteria extrema over all the time-steps - !> @{ - real(kind(0d0)) :: icfl_max !< ICFL criterion maximum - real(kind(0d0)) :: vcfl_max !< VCFL criterion maximum - real(kind(0d0)) :: ccfl_max !< CCFL criterion maximum - real(kind(0d0)) :: Rc_min !< Rc criterion maximum - !> @} procedure(s_write_abstract_data_files), pointer :: s_write_data_files => null() contains - !> The purpose of this subroutine is to open a new or pre- - !! existing run-time information file and append to it the - !! basic header information relevant to current simulation. - !! In general, this requires generating a table header for - !! those stability criteria which will be written at every - !! time-step. - subroutine s_open_run_time_information_file - - character(LEN=name_len) :: file_name = 'run_time.inf' !< - !! Name of the run-time information file - - character(LEN=path_len + name_len) :: file_path !< - !! Relative path to a file in the case directory - - character(LEN=8) :: file_date !< - !! Creation date of the run-time information file - - ! Opening the run-time information file - file_path = trim(case_dir)//'/'//trim(file_name) - - open (3, FILE=trim(file_path), & - FORM='formatted', & - STATUS='replace') - - write (3, '(A)') 'Description: Stability information at '// & - 'each time-step of the simulation. This' - write (3, '(13X,A)') 'data is composed of the inviscid '// & - 'Courant–Friedrichs–Lewy (ICFL)' - write (3, '(13X,A)') 'number, the viscous CFL (VCFL) number, '// & - 'the capillary CFL (CCFL)' - write (3, '(13X,A)') 'number and the cell Reynolds (Rc) '// & - 'number. Please note that only' - write (3, '(13X,A)') 'those stability conditions pertinent '// & - 'to the physics included in' - write (3, '(13X,A)') 'the current computation are displayed.' - - call date_and_time(DATE=file_date) - - write (3, '(A)') 'Date: '//file_date(5:6)//'/'// & - file_date(7:8)//'/'// & - file_date(3:4) - - write (3, '(A)') ''; write (3, '(A)') '' - - ! Generating table header for the stability criteria to be outputted - if (cfl_dt) then - if (any(Re_size > 0)) then - write (1, '(A)') '==== Time-steps ====== dt ===== Time ======= ICFL '// & - 'Max ==== VCFL Max ====== Rc Min =======' - else - write (1, '(A)') '=========== Time-steps ============== dt ===== Time '// & - '============== ICFL Max =============' - end if - else - if (any(Re_size > 0)) then - write (1, '(A)') '==== Time-steps ====== Time ======= ICFL '// & - 'Max ==== VCFL Max ====== Rc Min =======' - else - write (1, '(A)') '=========== Time-steps ============== Time '// & - '============== ICFL Max =============' - end if - end if - - end subroutine s_open_run_time_information_file - !> This opens a formatted data file where the root processor !! can write out flow probe information subroutine s_open_probe_files @@ -216,155 +119,6 @@ contains end subroutine s_open_probe_files - !> The goal of the procedure is to output to the run-time - !! information file the stability criteria extrema in the - !! entire computational domain and at the given time-step. - !! Moreover, the subroutine is also in charge of tracking - !! these stability criteria extrema over all time-steps. - !! @param q_prim_vf Cell-average primitive variables - !! @param t_step Current time step - subroutine s_write_run_time_information(q_prim_vf, t_step) - - type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf - integer, intent(IN) :: t_step - - real(kind(0d0)), dimension(num_fluids) :: alpha_rho !< Cell-avg. partial density - real(kind(0d0)) :: rho !< Cell-avg. density - real(kind(0d0)), dimension(num_dims) :: vel !< Cell-avg. velocity - real(kind(0d0)) :: vel_sum !< Cell-avg. velocity sum - real(kind(0d0)) :: pres !< Cell-avg. pressure - real(kind(0d0)), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction - real(kind(0d0)) :: gamma !< Cell-avg. sp. heat ratio - real(kind(0d0)) :: pi_inf !< Cell-avg. liquid stiffness function - real(kind(0d0)) :: qv !< Cell-avg. fluid reference energy - real(kind(0d0)) :: c !< Cell-avg. sound speed - real(kind(0d0)) :: E !< Cell-avg. energy - real(kind(0d0)) :: H !< Cell-avg. enthalpy - real(kind(0d0)), dimension(2) :: Re !< Cell-avg. Reynolds numbers - - ! ICFL, VCFL, CCFL and Rc stability criteria extrema for the current - ! time-step and located on both the local (loc) and the global (glb) - ! computational domains - - real(kind(0d0)) :: blkmod1, blkmod2 !< - !! Fluid bulk modulus for Woods mixture sound speed - - integer :: i, j, k, l, q !< Generic loop iterators - - integer :: Nfq - real(kind(0d0)) :: fltr_dtheta !< - !! Modified dtheta accounting for Fourier filtering in azimuthal direction. - - ! Computing Stability Criteria at Current Time-step ================ - !$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho, vel, alpha, Re, fltr_dtheta, Nfq) - do l = 0, p - do k = 0, n - do j = 0, m - - call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l) - - ! Compute mixture sound speed - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, c) - - if (any(Re_size > 0)) then - call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf, vcfl_sf, Rc_sf) - else - call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf) - end if - - end do - end do - end do - ! END: Computing Stability Criteria at Current Time-step =========== - - ! Determining local stability criteria extrema at current time-step - -#ifdef CRAY_ACC_WAR - !$acc update host(icfl_sf) - - if (any(Re_size > 0)) then - !$acc update host(vcfl_sf, Rc_sf) - end if - - icfl_max_loc = maxval(icfl_sf) - - if (any(Re_size > 0)) then - vcfl_max_loc = maxval(vcfl_sf) - Rc_min_loc = minval(Rc_sf) - end if -#else - !$acc kernels - icfl_max_loc = maxval(icfl_sf) - !$acc end kernels - - if (any(Re_size > 0)) then - !$acc kernels - vcfl_max_loc = maxval(vcfl_sf) - Rc_min_loc = minval(Rc_sf) - !$acc end kernels - end if -#endif - - ! Determining global stability criteria extrema at current time-step - if (num_procs > 1) then - call s_mpi_reduce_stability_criteria_extrema(icfl_max_loc, & - vcfl_max_loc, & - ccfl_max_loc, & - Rc_min_loc, & - icfl_max_glb, & - vcfl_max_glb, & - ccfl_max_glb, & - Rc_min_glb) - else - icfl_max_glb = icfl_max_loc - if (any(Re_size > 0)) vcfl_max_glb = vcfl_max_loc - if (any(Re_size > 0)) Rc_min_glb = Rc_min_loc - end if - - ! Determining the stability criteria extrema over all the time-steps - if (icfl_max_glb > icfl_max) icfl_max = icfl_max_glb - - if (any(Re_size > 0)) then - if (vcfl_max_glb > vcfl_max) vcfl_max = vcfl_max_glb - if (Rc_min_glb < Rc_min) Rc_min = Rc_min_glb - end if - - ! Outputting global stability criteria extrema at current time-step - if (proc_rank == 0) then - if (any(Re_size > 0)) then - write (1, '(6X,I8,F10.6,6X,6X,F10.6,6X,F9.6,6X,F9.6,6X,F10.6)') & - t_step, dt, t_step*dt, icfl_max_glb, & - vcfl_max_glb, & - Rc_min_glb - else - write (1, '(13X,I8,14X,F10.6,14X,F10.6,13X,F9.6)') & - t_step, dt, t_step*dt, icfl_max_glb - end if - - if (icfl_max_glb /= icfl_max_glb) then - call s_mpi_abort('ICFL is NaN. Exiting ...') - elseif (icfl_max_glb > 1d0) then - print *, 'icfl', icfl_max_glb - call s_mpi_abort('ICFL is greater than 1.0. Exiting ...') - end if - - do i = chemxb, chemxe - !@:ASSERT(all(q_prim_vf(i)%sf(:,:,:) >= -1d0), "bad conc") - !@:ASSERT(all(q_prim_vf(i)%sf(:,:,:) <= 2d0), "bad conc") - end do - - if (vcfl_max_glb /= vcfl_max_glb) then - call s_mpi_abort('VCFL is NaN. Exiting ...') - elseif (vcfl_max_glb > 1d0) then - print *, 'vcfl', vcfl_max_glb - call s_mpi_abort('VCFL is greater than 1.0. Exiting ...') - end if - end if - - call s_mpi_barrier() - - end subroutine s_write_run_time_information - !> The goal of this subroutine is to output the grid and !! conservative variables data files for given time-step. !! @param q_cons_vf Cell-average conservative variables @@ -1566,33 +1320,6 @@ contains end subroutine s_write_probe_files - !> The goal of this subroutine is to write to the run-time - !! information file basic footer information applicable to - !! the current computation and to close the file when done. - !! The footer contains the stability criteria extrema over - !! all of the time-steps and the simulation run-time. - subroutine s_close_run_time_information_file - - real(kind(0d0)) :: run_time !< Run-time of the simulation - ! Writing the footer of and closing the run-time information file - write (3, '(A)') '----------------------------------------'// & - '----------------------------------------' - write (3, '(A)') '' - - write (3, '(A,F9.6)') 'ICFL Max: ', icfl_max - if (any(Re_size > 0)) write (3, '(A,F9.6)') 'VCFL Max: ', vcfl_max - if (any(Re_size > 0)) write (3, '(A,F10.6)') 'Rc Min: ', Rc_min - - call cpu_time(run_time) - - write (3, '(A)') '' - write (3, '(A,I0,A)') 'Run-time: ', int(anint(run_time)), 's' - write (3, '(A)') '========================================'// & - '========================================' - close (3) - - end subroutine s_close_run_time_information_file - !> Closes probe files subroutine s_close_probe_files @@ -1613,18 +1340,6 @@ contains integer :: i !< Generic loop iterator - ! Allocating/initializing ICFL, VCFL, CCFL and Rc stability criteria - @:ALLOCATE_GLOBAL(icfl_sf(0:m, 0:n, 0:p)) - icfl_max = 0d0 - - if (any(Re_size > 0)) then - @:ALLOCATE_GLOBAL(vcfl_sf(0:m, 0:n, 0:p)) - @:ALLOCATE_GLOBAL(Rc_sf (0:m, 0:n, 0:p)) - - vcfl_max = 0d0 - Rc_min = 1d3 - end if - ! Associating the procedural pointer to the appropriate subroutine ! that will be utilized in the conversion to the mixture variables @@ -1652,12 +1367,6 @@ contains integer :: i !< Generic loop iterator - ! Deallocating the ICFL, VCFL, CCFL, and Rc stability criteria - @:DEALLOCATE_GLOBAL(icfl_sf) - if (any(Re_size > 0)) then - @:DEALLOCATE_GLOBAL(vcfl_sf, Rc_sf) - end if - ! Disassociating the pointer to the procedure that was utilized to ! to convert mixture or species variables to the mixture variables s_convert_to_mixture_variables => null() diff --git a/src/simulation/m_sim_helpers.f90 b/src/simulation/m_sim_helpers.f90 deleted file mode 100644 index be95197f6..000000000 --- a/src/simulation/m_sim_helpers.f90 +++ /dev/null @@ -1,250 +0,0 @@ -module m_sim_helpers - - ! Dependencies ============================================================= - use m_derived_types !< Definitions of the derived types - - use m_global_parameters - - use m_variables_conversion - - use m_mpi_proxy - ! ========================================================================== - - implicit none - - private; public :: s_compute_enthalpy, & - s_compute_stability_from_dt, & - s_compute_dt_from_cfl - -contains - - !> Computes enthalpy - !! @param q_prim_vf cell centered primitive variables - !! @param pres mixture pressure - !! @param rho mixture density - !! @param gamma mixture gamma - !! @param pi_inf mixture pi_inf - !! @param Re mixture reynolds number - !! @param H mixture enthalpy - !! @param alpha component alphas - !! @param vel directional velocities - !! @param vel_sum squard sum of velocity components - !! @param j x index - !! @param k y index - !! @param l z index - subroutine s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l) - !$acc routine seq - type(scalar_field), dimension(sys_size) :: q_prim_vf - real(kind(0d0)), dimension(num_fluids) :: alpha_rho - real(kind(0d0)), dimension(num_fluids) :: alpha - real(kind(0d0)), dimension(num_dims) :: vel - real(kind(0d0)) :: rho, gamma, pi_inf, qv, vel_sum, E, H, pres - real(kind(0d0)), dimension(2) :: Re - integer :: i, j, k, l - - do i = 1, num_fluids - alpha_rho(i) = q_prim_vf(i)%sf(j, k, l) - alpha(i) = q_prim_vf(E_idx + i)%sf(j, k, l) - end do - - if (bubbles) then - call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re, j, k, l) - else - call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re, j, k, l) - end if - - do i = 1, num_dims - vel(i) = q_prim_vf(contxe + i)%sf(j, k, l) - end do - - vel_sum = 0d0 - do i = 1, num_dims - vel_sum = vel_sum + vel(i)**2d0 - end do - - pres = q_prim_vf(E_idx)%sf(j, k, l) - - E = gamma*pres + pi_inf + 5d-1*rho*vel_sum + qv - - H = (E + pres)/rho - - end subroutine s_compute_enthalpy - - !> Computes stability criterion for a specified dt - !! @param vel directional velocities - !! @param c mixture speed of sound - !! @param Re_l mixture Reynolds number - !! @param j x index - !! @param k y index - !! @param l z index - !! @param icfl_sf cell centered inviscid cfl number - !! @param vcfl_sf (optional) cell centered viscous cfl number - !! @param Rc_sf (optional) cell centered Rc - subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl_sf, vcfl_sf, Rc_sf) - !$acc routine seq - real(kind(0d0)), dimension(num_dims) :: vel - real(kind(0d0)) :: c, icfl_dt, vcfl_dt, rho - real(kind(0d0)), dimension(0:m, 0:n, 0:p) :: icfl_sf - real(kind(0d0)), dimension(0:m, 0:n, 0:p), optional :: vcfl_sf, Rc_sf - real(kind(0d0)) :: fltr_dtheta !< - !! Modified dtheta accounting for Fourier filtering in azimuthal direction. - integer :: j, k, l - integer :: Nfq - real(kind(0d0)), dimension(2) :: Re_l - - if (grid_geometry == 3) then - if (k == 0) then - fltr_dtheta = 2d0*pi*y_cb(0)/3d0 - elseif (k <= fourier_rings) then - Nfq = min(floor(2d0*real(k, kind(0d0))*pi), (p + 1)/2 + 1) - fltr_dtheta = 2d0*pi*y_cb(k - 1)/real(Nfq, kind(0d0)) - else - fltr_dtheta = y_cb(k - 1)*dz(l) - end if - end if - - if (p > 0) then - !3D - if (grid_geometry == 3) then - icfl_sf(j, k, l) = dt/min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c), & - fltr_dtheta/(abs(vel(3)) + c)) - else - icfl_sf(j, k, l) = dt/min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c), & - dz(l)/(abs(vel(3)) + c)) - end if - - if (any(Re_size > 0)) then - - if (grid_geometry == 3) then - vcfl_sf(j, k, l) = maxval(dt/Re_l/rho) & - /min(dx(j), dy(k), fltr_dtheta)**2d0 - - Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), & - dy(k)*(abs(vel(2)) + c), & - fltr_dtheta*(abs(vel(3)) + c)) & - /maxval(1d0/Re_l) - else - vcfl_sf(j, k, l) = maxval(dt/Re_l/rho) & - /min(dx(j), dy(k), dz(l))**2d0 - - Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), & - dy(k)*(abs(vel(2)) + c), & - dz(l)*(abs(vel(3)) + c)) & - /maxval(1d0/Re_l) - end if - - end if - - elseif (n > 0) then - !2D - icfl_sf(j, k, l) = dt/min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c)) - - if (any(Re_size > 0)) then - - vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/min(dx(j), dy(k))**2d0 - - Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), & - dy(k)*(abs(vel(2)) + c)) & - /maxval(1d0/Re_l) - - end if - - else - !1D - icfl_sf(j, k, l) = (dt/dx(j))*(abs(vel(1)) + c) - - if (any(Re_size > 0)) then - - vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/dx(j)**2d0 - - Rc_sf(j, k, l) = dx(j)*(abs(vel(1)) + c)/maxval(1d0/Re_l) - - end if - - end if - - end subroutine s_compute_stability_from_dt - - !> Computes dt for a specified CFL number - !! @param vel directional velocities - !! @param max_dt cell centered maximum dt - !! @param rho cell centered density - !! @param Re_l cell centered Reynolds number - !! @param j x coordinate - !! @param k y coordinate - !! @param l z coordinate - subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l) - !$acc routine seq - real(kind(0d0)), dimension(num_dims) :: vel - real(kind(0d0)) :: c, icfl_dt, vcfl_dt, rho - real(kind(0d0)), dimension(0:m, 0:n, 0:p) :: max_dt - real(kind(0d0)) :: fltr_dtheta !< - !! Modified dtheta accounting for Fourier filtering in azimuthal direction. - integer :: j, k, l - integer :: Nfq - real(kind(0d0)), dimension(2) :: Re_l - - if (grid_geometry == 3) then - if (k == 0) then - fltr_dtheta = 2d0*pi*y_cb(0)/3d0 - elseif (k <= fourier_rings) then - Nfq = min(floor(2d0*real(k, kind(0d0))*pi), (p + 1)/2 + 1) - fltr_dtheta = 2d0*pi*y_cb(k - 1)/real(Nfq, kind(0d0)) - else - fltr_dtheta = y_cb(k - 1)*dz(l) - end if - end if - - if (p > 0) then - !3D - if (grid_geometry == 3) then - icfl_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c), & - fltr_dtheta/(abs(vel(3)) + c)) - else - icfl_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c), & - dz(l)/(abs(vel(3)) + c)) - end if - - if (any(Re_size > 0)) then - if (grid_geometry == 3) then - vcfl_dt = cfl_target*(min(dx(j), dy(k), fltr_dtheta)**2d0) & - /minval(1/(rho*Re_l)) - else - vcfl_dt = cfl_target*(min(dx(j), dy(k), dz(l))**2d0) & - /minval(1/(rho*Re_l)) - end if - end if - - elseif (n > 0) then - !2D - icfl_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c)) - - if (any(Re_size > 0)) then - vcfl_dt = cfl_target*(min(dx(j), dy(k))**2d0)/maxval((1/Re_l)/rho) - end if - - else - !1D - icfl_dt = cfl_target*(dx(j)/(abs(vel(1)) + c)) - - if (any(Re_size > 0)) then - vcfl_dt = cfl_target*(dx(j)**2d0)/minval(1/(rho*Re_l)) - end if - - end if - - if (any(re_size > 0)) then - max_dt(j, k, l) = min(icfl_dt, vcfl_dt) - else - max_dt(j, k, l) = icfl_dt - end if - - end subroutine s_compute_dt_from_cfl - -end module m_sim_helpers diff --git a/src/simulation/m_sim_helpers.fpp b/src/simulation/m_sim_helpers.fpp new file mode 100644 index 000000000..f533ff266 --- /dev/null +++ b/src/simulation/m_sim_helpers.fpp @@ -0,0 +1,1272 @@ +#:include 'macros.fpp' + +module m_sim_helpers + + ! Dependencies ============================================================= + use m_derived_types !< Definitions of the derived types + + use m_global_parameters + + use m_variables_conversion + + use m_mpi_proxy + + use m_data_output + + use m_compile_specific + ! ========================================================================== + + implicit none + + private; public :: s_compute_enthalpy, & + s_read_data_files, & + s_read_serial_data_files, & + s_read_parallel_data_files, & + s_compute_stability_from_dt, & + s_compute_dt_from_cfl, & + s_comprehensive_debug, & + s_open_run_time_information_file, & + s_write_run_time_information, & + s_close_run_time_information_file, & + s_initialize_sim_helpers_module, & + s_finalize_sim_helpers_module + + abstract interface ! =================================================== + + !! @param q_cons_vf Conservative variables + subroutine s_read_abstract_data_files(q_cons_vf, ib_markers) + + import :: scalar_field, integer_field, sys_size, pres_field + + type(scalar_field), & + dimension(sys_size), & + intent(inout) :: q_cons_vf + + type(integer_field) :: ib_markers + + end subroutine s_read_abstract_data_files + + end interface ! ======================================================== + + procedure(s_read_abstract_data_files), pointer :: s_read_data_files => null() + +#ifdef CRAY_ACC_WAR + @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), icfl_sf) + @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), vcfl_sf) + @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), ccfl_sf) + @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), Rc_sf) + !$acc declare link(icfl_sf, vcfl_sf, ccfl_sf, Rc_sf) +#else + real(kind(0d0)), allocatable, dimension(:, :, :) :: icfl_sf !< ICFL stability criterion + real(kind(0d0)), allocatable, dimension(:, :, :) :: vcfl_sf !< VCFL stability criterion + real(kind(0d0)), allocatable, dimension(:, :, :) :: ccfl_sf !< CCFL stability criterion + real(kind(0d0)), allocatable, dimension(:, :, :) :: Rc_sf !< Rc stability criterion + !$acc declare create(icfl_sf, vcfl_sf, ccfl_sf, Rc_sf) +#endif + + real(kind(0d0)) :: icfl_max_loc, icfl_max_glb !< ICFL stability extrema on local and global grids + real(kind(0d0)) :: vcfl_max_loc, vcfl_max_glb !< VCFL stability extrema on local and global grids + real(kind(0d0)) :: ccfl_max_loc, ccfl_max_glb !< CCFL stability extrema on local and global grids + real(kind(0d0)) :: Rc_min_loc, Rc_min_glb !< Rc stability extrema on local and global grids + !$acc declare create(icfl_max_loc, icfl_max_glb, vcfl_max_loc, vcfl_max_glb, ccfl_max_loc, ccfl_max_glb, Rc_min_loc, Rc_min_glb) + + !> @name ICFL, VCFL, CCFL and Rc stability criteria extrema over all the time-steps + !> @{ + real(kind(0d0)) :: icfl_max !< ICFL criterion maximum + real(kind(0d0)) :: vcfl_max !< VCFL criterion maximum + real(kind(0d0)) :: ccfl_max !< CCFL criterion maximum + real(kind(0d0)) :: Rc_min !< Rc criterion maximum + !> @} + +contains + + subroutine s_initialize_sim_helpers_module() + + ! Allocating/initializing ICFL, VCFL, CCFL and Rc stability criteria + @:ALLOCATE_GLOBAL(icfl_sf(0:m, 0:n, 0:p)) + icfl_max = 0d0 + + if (any(Re_size > 0)) then + @:ALLOCATE_GLOBAL(vcfl_sf(0:m, 0:n, 0:p)) + @:ALLOCATE_GLOBAL(Rc_sf (0:m, 0:n, 0:p)) + + vcfl_max = 0d0 + Rc_min = 1d3 + end if + + ! Associate pointers for serial or parallel I/O + if (parallel_io .neqv. .true.) then + s_read_data_files => s_read_serial_data_files + s_write_data_files => s_write_serial_data_files + else + s_read_data_files => s_read_parallel_data_files + s_write_data_files => s_write_parallel_data_files + end if + + end subroutine s_initialize_sim_helpers_module + + !> Computes enthalpy + !! @param q_prim_vf cell centered primitive variables + !! @param pres mixture pressure + !! @param rho mixture density + !! @param gamma mixture gamma + !! @param pi_inf mixture pi_inf + !! @param Re mixture reynolds number + !! @param H mixture enthalpy + !! @param alpha component alphas + !! @param vel directional velocities + !! @param vel_sum squard sum of velocity components + !! @param j x index + !! @param k y index + !! @param l z index + subroutine s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l) + !$acc routine seq + type(scalar_field), dimension(sys_size) :: q_prim_vf + real(kind(0d0)), dimension(num_fluids) :: alpha_rho + real(kind(0d0)), dimension(num_fluids) :: alpha + real(kind(0d0)), dimension(num_dims) :: vel + real(kind(0d0)) :: rho, gamma, pi_inf, qv, vel_sum, E, H, pres + real(kind(0d0)), dimension(2) :: Re + integer :: i, j, k, l + + do i = 1, num_fluids + alpha_rho(i) = q_prim_vf(i)%sf(j, k, l) + alpha(i) = q_prim_vf(E_idx + i)%sf(j, k, l) + end do + + if (bubbles) then + call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re, j, k, l) + else + call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re, j, k, l) + end if + + do i = 1, num_dims + vel(i) = q_prim_vf(contxe + i)%sf(j, k, l) + end do + + vel_sum = 0d0 + do i = 1, num_dims + vel_sum = vel_sum + vel(i)**2d0 + end do + + pres = q_prim_vf(E_idx)%sf(j, k, l) + + E = gamma*pres + pi_inf + 5d-1*rho*vel_sum + qv + + H = (E + pres)/rho + + end subroutine s_compute_enthalpy + + !> Computes stability criterion for a specified dt + !! @param vel directional velocities + !! @param c mixture speed of sound + !! @param Re_l mixture Reynolds number + !! @param j x index + !! @param k y index + !! @param l z index + !! @param icfl_sf cell centered inviscid cfl number + !! @param vcfl_sf (optional) cell centered viscous cfl number + !! @param Rc_sf (optional) cell centered Rc + subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl_sf, vcfl_sf, Rc_sf) + !$acc routine seq + real(kind(0d0)), dimension(num_dims) :: vel + real(kind(0d0)) :: c, icfl_dt, vcfl_dt, rho + real(kind(0d0)), dimension(0:m, 0:n, 0:p) :: icfl_sf + real(kind(0d0)), dimension(0:m, 0:n, 0:p), optional :: vcfl_sf, Rc_sf + real(kind(0d0)) :: fltr_dtheta !< + !! Modified dtheta accounting for Fourier filtering in azimuthal direction. + integer :: j, k, l + integer :: Nfq + real(kind(0d0)), dimension(2) :: Re_l + + if (grid_geometry == 3) then + if (k == 0) then + fltr_dtheta = 2d0*pi*y_cb(0)/3d0 + elseif (k <= fourier_rings) then + Nfq = min(floor(2d0*real(k, kind(0d0))*pi), (p + 1)/2 + 1) + fltr_dtheta = 2d0*pi*y_cb(k - 1)/real(Nfq, kind(0d0)) + else + fltr_dtheta = y_cb(k - 1)*dz(l) + end if + end if + + if (p > 0) then + !3D + if (grid_geometry == 3) then + icfl_sf(j, k, l) = dt/min(dx(j)/(abs(vel(1)) + c), & + dy(k)/(abs(vel(2)) + c), & + fltr_dtheta/(abs(vel(3)) + c)) + else + icfl_sf(j, k, l) = dt/min(dx(j)/(abs(vel(1)) + c), & + dy(k)/(abs(vel(2)) + c), & + dz(l)/(abs(vel(3)) + c)) + end if + + if (any(Re_size > 0)) then + + if (grid_geometry == 3) then + vcfl_sf(j, k, l) = maxval(dt/Re_l/rho) & + /min(dx(j), dy(k), fltr_dtheta)**2d0 + + Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), & + dy(k)*(abs(vel(2)) + c), & + fltr_dtheta*(abs(vel(3)) + c)) & + /maxval(1d0/Re_l) + else + vcfl_sf(j, k, l) = maxval(dt/Re_l/rho) & + /min(dx(j), dy(k), dz(l))**2d0 + + Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), & + dy(k)*(abs(vel(2)) + c), & + dz(l)*(abs(vel(3)) + c)) & + /maxval(1d0/Re_l) + end if + + end if + + elseif (n > 0) then + !2D + icfl_sf(j, k, l) = dt/min(dx(j)/(abs(vel(1)) + c), & + dy(k)/(abs(vel(2)) + c)) + + if (any(Re_size > 0)) then + + vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/min(dx(j), dy(k))**2d0 + + Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), & + dy(k)*(abs(vel(2)) + c)) & + /maxval(1d0/Re_l) + + end if + + else + !1D + icfl_sf(j, k, l) = (dt/dx(j))*(abs(vel(1)) + c) + + if (any(Re_size > 0)) then + + vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/dx(j)**2d0 + + Rc_sf(j, k, l) = dx(j)*(abs(vel(1)) + c)/maxval(1d0/Re_l) + + end if + + end if + + end subroutine s_compute_stability_from_dt + + !> Computes dt for a specified CFL number + !! @param vel directional velocities + !! @param max_dt cell centered maximum dt + !! @param rho cell centered density + !! @param Re_l cell centered Reynolds number + !! @param j x coordinate + !! @param k y coordinate + !! @param l z coordinate + subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l) + !$acc routine seq + real(kind(0d0)), dimension(num_dims) :: vel + real(kind(0d0)) :: c, icfl_dt, vcfl_dt, rho + real(kind(0d0)), dimension(0:m, 0:n, 0:p) :: max_dt + real(kind(0d0)) :: fltr_dtheta !< + !! Modified dtheta accounting for Fourier filtering in azimuthal direction. + integer :: j, k, l + integer :: Nfq + real(kind(0d0)), dimension(2) :: Re_l + + if (grid_geometry == 3) then + if (k == 0) then + fltr_dtheta = 2d0*pi*y_cb(0)/3d0 + elseif (k <= fourier_rings) then + Nfq = min(floor(2d0*real(k, kind(0d0))*pi), (p + 1)/2 + 1) + fltr_dtheta = 2d0*pi*y_cb(k - 1)/real(Nfq, kind(0d0)) + else + fltr_dtheta = y_cb(k - 1)*dz(l) + end if + end if + + if (p > 0) then + !3D + if (grid_geometry == 3) then + icfl_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), & + dy(k)/(abs(vel(2)) + c), & + fltr_dtheta/(abs(vel(3)) + c)) + else + icfl_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), & + dy(k)/(abs(vel(2)) + c), & + dz(l)/(abs(vel(3)) + c)) + end if + + if (any(Re_size > 0)) then + if (grid_geometry == 3) then + vcfl_dt = cfl_target*(min(dx(j), dy(k), fltr_dtheta)**2d0) & + /minval(1/(rho*Re_l)) + else + vcfl_dt = cfl_target*(min(dx(j), dy(k), dz(l))**2d0) & + /minval(1/(rho*Re_l)) + end if + end if + + elseif (n > 0) then + !2D + icfl_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), & + dy(k)/(abs(vel(2)) + c)) + + if (any(Re_size > 0)) then + vcfl_dt = cfl_target*(min(dx(j), dy(k))**2d0)/maxval((1/Re_l)/rho) + end if + + else + !1D + icfl_dt = cfl_target*(dx(j)/(abs(vel(1)) + c)) + + if (any(Re_size > 0)) then + vcfl_dt = cfl_target*(dx(j)**2d0)/minval(1/(rho*Re_l)) + end if + + end if + + if (any(re_size > 0)) then + max_dt(j, k, l) = min(icfl_dt, vcfl_dt) + else + max_dt(j, k, l) = icfl_dt + end if + + end subroutine s_compute_dt_from_cfl + + subroutine s_comprehensive_debug(q_cons_vf, q_prim_vf, t_step, stage) + + type(scalar_field), dimension(sys_size) :: q_cons_vf, q_prim_vf + integer, intent(in) :: t_step, stage + integer :: j, k, l, i + integer errors + logical :: exists + + character(LEN=name_len) :: file_name = 'comp_debug.txt' + character(LEN=path_len + name_len) :: file_path + character(100) :: str_format + + ! Opening the run-time information file + file_path = trim(case_dir)//'/'//trim(file_name) + + str_format = "(I9, A, I3, A, I4, I4, I4, A, I2, A, I5, A, I5, I5, I5)" + + open (12, FILE=trim(file_path), & + STATUS='replace') + + errors = 0 + + ! Check all variables for NaNs + do i = 1, sys_size + !$acc update host(q_cons_vf(i)%sf) + do l = 0, p + do k = 0, n + do j = 0, m + if (ieee_is_nan(q_cons_vf(i)%sf(j, k, l))) then + write(12, str_format) t_step, " NaN(s) in conservative variables after RK stage ", & + stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & + " (m, n, p)", m, n, p + errors = errors + 1 + end if + end do + end do + end do + end do + + ! Check for invalid volume fractions + do i = advxb, advxe + do l = 0, p + do k = 0, n + do j = 0, m + if (q_cons_vf(i)%sf(j, k, l) < 0d0) then + write(12, str_format) t_step, " Volume fraction < 0 after RK stage ", & + stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & + " (m, n, p)", m, n, p + errors = errors + 1 + elseif (q_cons_vf(i)%sf(j, k, l) > 1d0 + verysmall) then + write(12, str_format) t_step, " Volume fraction > 1 after RK stage ", & + stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & + " (m, n, p)", m, n, p + errors = errors + 1 + end if + end do + end do + end do + end do + + ! Check for invalid densities + do i = contxb, contxe + do l = 0, p + do k = 0, n + do j = 0, m + if (q_cons_vf(advxb + i -1)%sf(j, k, l) < 0d0 .and. q_cons_vf(i)%sf(j, k, l) < 0d0 .or. & + q_cons_vf(advxb + i -1)%sf(j, k, l) > 0d0 .and. q_cons_Vf(i)%sf(j, k, l) < 0d0) then + print*, q_cons_vf(advxb + i - 1)%sf(j, k, l), q_cons_vf(i)%sf(j, k, l) + write(12, str_format) t_step, " Density is negative after RK stage ", & + stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & + " (m, n, p)", m, n, p + errors = errors + 1 + end if + end do + end do + end do + end do + + if (errors /= 0) then + close(12) + call s_write_data_files(q_cons_vf, q_prim_vf, t_step) + call s_mpi_abort("Errors found in conservative variables") + endif + + write(12, "(I3)") -1 + close(12) + + end subroutine s_comprehensive_debug + + !! initial condition and grid data files. The cell-average + !! conservative variables constitute the former, while the + !! cell-boundary locations in x-, y- and z-directions make + !! up the latter. This procedure also calculates the cell- + !! width distributions from the cell-boundary locations. + !! @param q_cons_vf Cell-averaged conservative variables + subroutine s_read_serial_data_files(q_cons_vf, ib_markers) + + type(scalar_field), dimension(sys_size), intent(INOUT) :: q_cons_vf + type(integer_field) :: ib_markers + + character(LEN=path_len + 2*name_len) :: t_step_dir !< + !! Relative path to the starting time-step directory + + character(LEN=path_len + 3*name_len) :: file_path !< + !! Relative path to the grid and conservative variables data files + + logical :: file_exist !< + ! Logical used to check the existence of the data files + + integer :: i, r !< Generic loop iterator + + ! Confirming that the directory from which the initial condition and + ! the grid data files are to be read in exists and exiting otherwise + if (cfl_dt) then + write (t_step_dir, '(A,I0,A,I0)') & + trim(case_dir)//'/p_all/p', proc_rank, '/', n_start + else + write (t_step_dir, '(A,I0,A,I0)') & + trim(case_dir)//'/p_all/p', proc_rank, '/', t_step_start + end if + + file_path = trim(t_step_dir)//'/.' + call my_inquire(file_path, file_exist) + + if (file_exist .neqv. .true.) then + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + + ! Cell-boundary Locations in x-direction =========================== + file_path = trim(t_step_dir)//'/x_cb.dat' + + inquire (FILE=trim(file_path), EXIST=file_exist) + + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) x_cb(-1:m); close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + + dx(0:m) = x_cb(0:m) - x_cb(-1:m - 1) + x_cc(0:m) = x_cb(-1:m - 1) + dx(0:m)/2d0 + + if (ib) then + do i = 1, num_ibs + if (patch_ib(i)%c > 0) then + Np = int((patch_ib(i)%p*patch_ib(i)%c/dx(0))*20) + int(((patch_ib(i)%c - patch_ib(i)%p*patch_ib(i)%c)/dx(0))*20) + 1 + end if + end do + end if + ! ================================================================== + + ! Cell-boundary Locations in y-direction =========================== + if (n > 0) then + + file_path = trim(t_step_dir)//'/y_cb.dat' + + inquire (FILE=trim(file_path), EXIST=file_exist) + + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) y_cb(-1:n); close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + + dy(0:n) = y_cb(0:n) - y_cb(-1:n - 1) + y_cc(0:n) = y_cb(-1:n - 1) + dy(0:n)/2d0 + + end if + ! ================================================================== + + ! Cell-boundary Locations in z-direction =========================== + if (p > 0) then + + file_path = trim(t_step_dir)//'/z_cb.dat' + + inquire (FILE=trim(file_path), EXIST=file_exist) + + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) z_cb(-1:p); close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + + dz(0:p) = z_cb(0:p) - z_cb(-1:p - 1) + z_cc(0:p) = z_cb(-1:p - 1) + dz(0:p)/2d0 + + end if + ! ================================================================== + + do i = 1, sys_size + write (file_path, '(A,I0,A)') & + trim(t_step_dir)//'/q_cons_vf', i, '.dat' + inquire (FILE=trim(file_path), EXIST=file_exist) + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) q_cons_vf(i)%sf(0:m, 0:n, 0:p); close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + end do + + if ((bubbles .eqv. .true.) .or. (hypoelasticity .eqv. .true.)) then + ! Read pb and mv for non-polytropic qbmm + if (qbmm .and. .not. polytropic) then + do i = 1, nb + do r = 1, nnode + write (file_path, '(A,I0,A)') & + trim(t_step_dir)//'/pb', sys_size + (i - 1)*nnode + r, '.dat' + inquire (FILE=trim(file_path), EXIST=file_exist) + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) pb_ts(1)%sf(0:m, 0:n, 0:p, r, i); close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + end do + end do + do i = 1, nb + do r = 1, nnode + write (file_path, '(A,I0,A)') & + trim(t_step_dir)//'/mv', sys_size + (i - 1)*nnode + r, '.dat' + inquire (FILE=trim(file_path), EXIST=file_exist) + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) mv_ts(1)%sf(0:m, 0:n, 0:p, r, i); close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + end do + end do + end if + end if + ! ================================================================== + + ! Read IBM Data ==================================================== + + if (ib) then + do i = 1, num_ibs + write (file_path, '(A,I0,A)') & + trim(t_step_dir)//'/ib.dat' + inquire (FILE=trim(file_path), EXIST=file_exist) + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) ib_markers%sf(0:m, 0:n, 0:p); close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + + if (patch_ib(i)%c > 0) then + + print *, "HERE Np", Np + allocate (airfoil_grid_u(1:Np)) + allocate (airfoil_grid_l(1:Np)) + + write (file_path, '(A)') & + trim(t_step_dir)//'/airfoil_u.dat' + inquire (FILE=trim(file_path), EXIST=file_exist) + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) airfoil_grid_u; close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + + write (file_path, '(A)') & + trim(t_step_dir)//'/airfoil_l.dat' + inquire (FILE=trim(file_path), EXIST=file_exist) + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) airfoil_grid_l; close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + end if + end do + + end if + + end subroutine s_read_serial_data_files + + !! @param q_cons_vf Conservative variables + subroutine s_read_parallel_data_files(q_cons_vf, ib_markers) + + type(scalar_field), & + dimension(sys_size), & + intent(INOUT) :: q_cons_vf + type(integer_field) :: ib_markers + +#ifdef MFC_MPI + + real(kind(0d0)), allocatable, dimension(:) :: x_cb_glb, y_cb_glb, z_cb_glb + + 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) :: 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_start_string + + integer :: i, j + + allocate (x_cb_glb(-1:m_glb)) + allocate (y_cb_glb(-1:n_glb)) + allocate (z_cb_glb(-1:p_glb)) + + ! Read in cell boundary locations in x-direction + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//'x_cb.dat' + inquire (FILE=trim(file_loc), EXIST=file_exist) + + if (file_exist) then + data_size = m_glb + 2 + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + call MPI_FILE_READ(ifile, x_cb_glb, data_size, MPI_DOUBLE_PRECISION, status, ierr) + call MPI_FILE_CLOSE(ifile, ierr) + else + call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') + end if + + ! Assigning local cell boundary locations + x_cb(-1:m) = x_cb_glb((start_idx(1) - 1):(start_idx(1) + m)) + ! Computing the cell width distribution + dx(0:m) = x_cb(0:m) - x_cb(-1:m - 1) + ! Computing the cell center locations + x_cc(0:m) = x_cb(-1:m - 1) + dx(0:m)/2d0 + + if (ib) then + do i = 1, num_ibs + if (patch_ib(i)%c > 0) then + Np = int((patch_ib(i)%p*patch_ib(i)%c/dx(0))*20) + int(((patch_ib(i)%c - patch_ib(i)%p*patch_ib(i)%c)/dx(0))*20) + 1 + allocate (MPI_IO_airfoil_IB_DATA%var(1:2*Np)) + print *, "HERE Np", Np + end if + end do + end if + + if (n > 0) then + ! Read in cell boundary locations in y-direction + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//'y_cb.dat' + inquire (FILE=trim(file_loc), EXIST=file_exist) + + if (file_exist) then + data_size = n_glb + 2 + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + call MPI_FILE_READ(ifile, y_cb_glb, data_size, MPI_DOUBLE_PRECISION, status, ierr) + call MPI_FILE_CLOSE(ifile, ierr) + else + call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') + end if + + ! Assigning local cell boundary locations + y_cb(-1:n) = y_cb_glb((start_idx(2) - 1):(start_idx(2) + n)) + ! Computing the cell width distribution + dy(0:n) = y_cb(0:n) - y_cb(-1:n - 1) + ! Computing the cell center locations + y_cc(0:n) = y_cb(-1:n - 1) + dy(0:n)/2d0 + + if (p > 0) then + ! Read in cell boundary locations in z-direction + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//'z_cb.dat' + inquire (FILE=trim(file_loc), EXIST=file_exist) + + if (file_exist) then + data_size = p_glb + 2 + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + call MPI_FILE_READ(ifile, z_cb_glb, data_size, MPI_DOUBLE_PRECISION, status, ierr) + call MPI_FILE_CLOSE(ifile, ierr) + else + call s_mpi_abort('File '//trim(file_loc)//'is missing. Exiting...') + end if + + ! Assigning local cell boundary locations + z_cb(-1:p) = z_cb_glb((start_idx(3) - 1):(start_idx(3) + p)) + ! Computing the cell width distribution + dz(0:p) = z_cb(0:p) - z_cb(-1:p - 1) + ! Computing the cell center locations + z_cc(0:p) = z_cb(-1:p - 1) + dz(0:p)/2d0 + + end if + end if + + if (file_per_process) then + if (cfl_dt) then + call s_int_to_str(n_start, t_step_start_string) + write (file_loc, '(I0,A1,I7.7,A)') n_start, '_', proc_rank, '.dat' + else + call s_int_to_str(t_step_start, t_step_start_string) + write (file_loc, '(I0,A1,I7.7,A)') t_step_start, '_', proc_rank, '.dat' + end if + + file_loc = trim(case_dir)//'/restart_data/lustre_'//trim(t_step_start_string)//trim(mpiiofs)//trim(file_loc) + inquire (FILE=trim(file_loc), EXIST=file_exist) + + if (file_exist) then + call MPI_FILE_OPEN(MPI_COMM_SELF, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + + ! Initialize MPI data I/O + + if (ib) then + call s_initialize_mpi_data(q_cons_vf, ib_markers) + else + call s_initialize_mpi_data(q_cons_vf) + end if + + ! Size of local arrays + data_size = (m + 1)*(n + 1)*(p + 1) + + ! Resize some integers so MPI can read even the biggest file + m_MOK = int(m_glb + 1, MPI_OFFSET_KIND) + n_MOK = int(n_glb + 1, MPI_OFFSET_KIND) + p_MOK = int(p_glb + 1, MPI_OFFSET_KIND) + WP_MOK = int(8d0, MPI_OFFSET_KIND) + MOK = int(1d0, MPI_OFFSET_KIND) + str_MOK = int(name_len, MPI_OFFSET_KIND) + NVARS_MOK = int(sys_size, MPI_OFFSET_KIND) + + ! Read the data for each variable + if (bubbles .or. hypoelasticity) then + + do i = 1, sys_size!adv_idx%end + var_MOK = int(i, MPI_OFFSET_KIND) + + call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & + MPI_DOUBLE_PRECISION, status, ierr) + end do + !Read pb and mv for non-polytropic qbmm + if (qbmm .and. .not. polytropic) then + do i = sys_size + 1, sys_size + 2*nb*nnode + var_MOK = int(i, MPI_OFFSET_KIND) + + call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & + MPI_DOUBLE_PRECISION, status, ierr) + end do + end if + else + do i = 1, adv_idx%end + var_MOK = int(i, MPI_OFFSET_KIND) + + call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & + MPI_DOUBLE_PRECISION, status, ierr) + end do + end if + + call s_mpi_barrier() + + call MPI_FILE_CLOSE(ifile, ierr) + + if (ib) then + + write (file_loc, '(A)') 'ib.dat' + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) + inquire (FILE=trim(file_loc), EXIST=file_exist) + + if (file_exist) then + + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + + disp = 0 + + call MPI_FILE_SET_VIEW(ifile, disp, MPI_INTEGER, MPI_IO_IB_DATA%view, & + 'native', mpi_info_int, ierr) + call MPI_FILE_READ(ifile, MPI_IO_IB_DATA%var%sf, data_size, & + MPI_INTEGER, status, ierr) + + else + call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') + end if + + end if + + else + call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') + end if + else + ! Open the file to read conservative variables + if (cfl_dt) then + write (file_loc, '(I0,A)') n_start, '.dat' + else + write (file_loc, '(I0,A)') t_step_start, '.dat' + end if + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) + inquire (FILE=trim(file_loc), EXIST=file_exist) + + if (file_exist) then + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + + ! Initialize MPI data I/O + + if (ib) then + call s_initialize_mpi_data(q_cons_vf, ib_markers) + else + + call s_initialize_mpi_data(q_cons_vf) + + end if + + + ! Size of local arrays + data_size = (m + 1)*(n + 1)*(p + 1) + + ! Resize some integers so MPI can read even the biggest file + m_MOK = int(m_glb + 1, MPI_OFFSET_KIND) + n_MOK = int(n_glb + 1, MPI_OFFSET_KIND) + p_MOK = int(p_glb + 1, MPI_OFFSET_KIND) + WP_MOK = int(8d0, MPI_OFFSET_KIND) + MOK = int(1d0, MPI_OFFSET_KIND) + str_MOK = int(name_len, MPI_OFFSET_KIND) + NVARS_MOK = int(sys_size, MPI_OFFSET_KIND) + + ! Read the data for each variable + if (bubbles .or. hypoelasticity) then + + do i = 1, sys_size!adv_idx%end + var_MOK = int(i, MPI_OFFSET_KIND) + ! Initial displacement to skip at beginning of file + disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1) + + call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_DATA%view(i), & + 'native', mpi_info_int, ierr) + call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & + MPI_DOUBLE_PRECISION, status, ierr) + end do + !Read pb and mv for non-polytropic qbmm + if (qbmm .and. .not. polytropic) then + do i = sys_size + 1, sys_size + 2*nb*nnode + var_MOK = int(i, MPI_OFFSET_KIND) + ! Initial displacement to skip at beginning of file + disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1) + + call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_DATA%view(i), & + 'native', mpi_info_int, ierr) + call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & + MPI_DOUBLE_PRECISION, status, ierr) + end do + end if + else + do i = 1, sys_size + var_MOK = int(i, MPI_OFFSET_KIND) + + ! Initial displacement to skip at beginning of file + disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1) + + call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_DATA%view(i), & + 'native', mpi_info_int, ierr) + call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & + MPI_DOUBLE_PRECISION, status, ierr) + + end do + end if + + call s_mpi_barrier() + + call MPI_FILE_CLOSE(ifile, ierr) + + if (ib) then + + write (file_loc, '(A)') 'ib.dat' + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) + inquire (FILE=trim(file_loc), EXIST=file_exist) + + if (file_exist) then + + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + + disp = 0 + + call MPI_FILE_SET_VIEW(ifile, disp, MPI_INTEGER, MPI_IO_IB_DATA%view, & + 'native', mpi_info_int, ierr) + call MPI_FILE_READ(ifile, MPI_IO_IB_DATA%var%sf, data_size, & + MPI_INTEGER, status, ierr) + + else + call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') + end if + + end if + + else + call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') + end if + + end if + + if (ib) then + + do j = 1, num_ibs + if (patch_ib(j)%c > 0) then + + print *, "HERE Np", Np + + allocate (airfoil_grid_u(1:Np)) + allocate (airfoil_grid_l(1:Np)) + + write (file_loc, '(A)') 'airfoil_l.dat' + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) + inquire (FILE=trim(file_loc), EXIST=file_exist) + if (file_exist) then + + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + + ! Initial displacement to skip at beginning of file + disp = 0 + + call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_airfoil_IB_DATA%view(1), & + 'native', mpi_info_int, ierr) + call MPI_FILE_READ(ifile, MPI_IO_airfoil_IB_DATA%var(1:Np), 3*Np, & + MPI_DOUBLE_PRECISION, status, ierr) + + end if + + write (file_loc, '(A)') 'airfoil_u.dat' + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) + inquire (FILE=trim(file_loc), EXIST=file_exist) + if (file_exist) then + + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + + ! Initial displacement to skip at beginning of file + disp = 0 + + call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_airfoil_IB_DATA%view(2), & + 'native', mpi_info_int, ierr) + call MPI_FILE_READ(ifile, MPI_IO_airfoil_IB_DATA%var(Np + 1:2*Np), 3*Np, & + MPI_DOUBLE_PRECISION, status, ierr) + end if + + do i = 1, Np + airfoil_grid_l(i)%x = MPI_IO_airfoil_IB_DATA%var(i)%x + airfoil_grid_l(i)%y = MPI_IO_airfoil_IB_DATA%var(i)%y + end do + + do i = 1, Np + airfoil_grid_u(i)%x = MPI_IO_airfoil_IB_DATA%var(Np + i)%x + airfoil_grid_u(i)%y = MPI_IO_airfoil_IB_DATA%var(Np + i)%y + end do + + end if + end do + end if + + deallocate (x_cb_glb, y_cb_glb, z_cb_glb) + +#endif + + end subroutine s_read_parallel_data_files + + !> The purpose of this subroutine is to open a new or pre- + !! existing run-time information file and append to it the + !! basic header information relevant to current simulation. + !! In general, this requires generating a table header for + !! those stability criteria which will be written at every + !! time-step. + subroutine s_open_run_time_information_file + + character(LEN=name_len) :: file_name = 'run_time.inf' !< + !! Name of the run-time information file + + character(LEN=path_len + name_len) :: file_path !< + !! Relative path to a file in the case directory + + character(LEN=8) :: file_date !< + !! Creation date of the run-time information file + + ! Opening the run-time information file + file_path = trim(case_dir)//'/'//trim(file_name) + + open (3, FILE=trim(file_path), & + FORM='formatted', & + STATUS='replace') + + write (3, '(A)') 'Description: Stability information at '// & + 'each time-step of the simulation. This' + write (3, '(13X,A)') 'data is composed of the inviscid '// & + 'Courant–Friedrichs–Lewy (ICFL)' + write (3, '(13X,A)') 'number, the viscous CFL (VCFL) number, '// & + 'the capillary CFL (CCFL)' + write (3, '(13X,A)') 'number and the cell Reynolds (Rc) '// & + 'number. Please note that only' + write (3, '(13X,A)') 'those stability conditions pertinent '// & + 'to the physics included in' + write (3, '(13X,A)') 'the current computation are displayed.' + + call date_and_time(DATE=file_date) + + write (3, '(A)') 'Date: '//file_date(5:6)//'/'// & + file_date(7:8)//'/'// & + file_date(3:4) + + write (3, '(A)') ''; write (3, '(A)') '' + + ! Generating table header for the stability criteria to be outputted + if (cfl_dt) then + if (any(Re_size > 0)) then + write (1, '(A)') '==== Time-steps ====== dt ===== Time ======= ICFL '// & + 'Max ==== VCFL Max ====== Rc Min =======' + else + write (1, '(A)') '=========== Time-steps ============== dt ===== Time '// & + '============== ICFL Max =============' + end if + else + if (any(Re_size > 0)) then + write (1, '(A)') '==== Time-steps ====== Time ======= ICFL '// & + 'Max ==== VCFL Max ====== Rc Min =======' + else + write (1, '(A)') '=========== Time-steps ============== Time '// & + '============== ICFL Max =============' + end if + end if + + end subroutine s_open_run_time_information_file + + !> The goal of the procedure is to output to the run-time + !! information file the stability criteria extrema in the + !! entire computational domain and at the given time-step. + !! Moreover, the subroutine is also in charge of tracking + !! these stability criteria extrema over all time-steps. + !! @param q_prim_vf Cell-average primitive variables + !! @param t_step Current time step + subroutine s_write_run_time_information(q_prim_vf, t_step) + + type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf + integer, intent(IN) :: t_step + + real(kind(0d0)), dimension(num_fluids) :: alpha_rho !< Cell-avg. partial density + real(kind(0d0)) :: rho !< Cell-avg. density + real(kind(0d0)), dimension(num_dims) :: vel !< Cell-avg. velocity + real(kind(0d0)) :: vel_sum !< Cell-avg. velocity sum + real(kind(0d0)) :: pres !< Cell-avg. pressure + real(kind(0d0)), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction + real(kind(0d0)) :: gamma !< Cell-avg. sp. heat ratio + real(kind(0d0)) :: pi_inf !< Cell-avg. liquid stiffness function + real(kind(0d0)) :: qv !< Cell-avg. fluid reference energy + real(kind(0d0)) :: c !< Cell-avg. sound speed + real(kind(0d0)) :: E !< Cell-avg. energy + real(kind(0d0)) :: H !< Cell-avg. enthalpy + real(kind(0d0)), dimension(2) :: Re !< Cell-avg. Reynolds numbers + + ! ICFL, VCFL, CCFL and Rc stability criteria extrema for the current + ! time-step and located on both the local (loc) and the global (glb) + ! computational domains + + real(kind(0d0)) :: blkmod1, blkmod2 !< + !! Fluid bulk modulus for Woods mixture sound speed + + integer :: i, j, k, l, q !< Generic loop iterators + + integer :: Nfq + real(kind(0d0)) :: fltr_dtheta !< + !! Modified dtheta accounting for Fourier filtering in azimuthal direction. + + ! Computing Stability Criteria at Current Time-step ================ + !$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho, vel, alpha, Re, fltr_dtheta, Nfq) + do l = 0, p + do k = 0, n + do j = 0, m + + call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l) + + ! Compute mixture sound speed + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, c) + + if (any(Re_size > 0)) then + call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf, vcfl_sf, Rc_sf) + else + call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf) + end if + + end do + end do + end do + ! END: Computing Stability Criteria at Current Time-step =========== + + ! Determining local stability criteria extrema at current time-step + +#ifdef CRAY_ACC_WAR + !$acc update host(icfl_sf) + + if (any(Re_size > 0)) then + !$acc update host(vcfl_sf, Rc_sf) + end if + + icfl_max_loc = maxval(icfl_sf) + + if (any(Re_size > 0)) then + vcfl_max_loc = maxval(vcfl_sf) + Rc_min_loc = minval(Rc_sf) + end if +#else + !$acc kernels + icfl_max_loc = maxval(icfl_sf) + !$acc end kernels + + if (any(Re_size > 0)) then + !$acc kernels + vcfl_max_loc = maxval(vcfl_sf) + Rc_min_loc = minval(Rc_sf) + !$acc end kernels + end if +#endif + + ! Determining global stability criteria extrema at current time-step + if (num_procs > 1) then + call s_mpi_reduce_stability_criteria_extrema(icfl_max_loc, & + vcfl_max_loc, & + ccfl_max_loc, & + Rc_min_loc, & + icfl_max_glb, & + vcfl_max_glb, & + ccfl_max_glb, & + Rc_min_glb) + else + icfl_max_glb = icfl_max_loc + if (any(Re_size > 0)) vcfl_max_glb = vcfl_max_loc + if (any(Re_size > 0)) Rc_min_glb = Rc_min_loc + end if + + ! Determining the stability criteria extrema over all the time-steps + if (icfl_max_glb > icfl_max) icfl_max = icfl_max_glb + + if (any(Re_size > 0)) then + if (vcfl_max_glb > vcfl_max) vcfl_max = vcfl_max_glb + if (Rc_min_glb < Rc_min) Rc_min = Rc_min_glb + end if + + ! Outputting global stability criteria extrema at current time-step + if (proc_rank == 0) then + if (any(Re_size > 0)) then + write (1, '(6X,I8,F10.6,6X,6X,F10.6,6X,F9.6,6X,F9.6,6X,F10.6)') & + t_step, dt, t_step*dt, icfl_max_glb, & + vcfl_max_glb, & + Rc_min_glb + else + write (1, '(13X,I8,14X,F10.6,14X,F10.6,13X,F9.6)') & + t_step, dt, t_step*dt, icfl_max_glb + end if + + if (icfl_max_glb /= icfl_max_glb) then + call s_mpi_abort('ICFL is NaN. Exiting ...') + elseif (icfl_max_glb > 1d0) then + print *, 'icfl', icfl_max_glb + call s_mpi_abort('ICFL is greater than 1.0. Exiting ...') + end if + + do i = chemxb, chemxe + !@:ASSERT(all(q_prim_vf(i)%sf(:,:,:) >= -1d0), "bad conc") + !@:ASSERT(all(q_prim_vf(i)%sf(:,:,:) <= 2d0), "bad conc") + end do + + if (vcfl_max_glb /= vcfl_max_glb) then + call s_mpi_abort('VCFL is NaN. Exiting ...') + elseif (vcfl_max_glb > 1d0) then + print *, 'vcfl', vcfl_max_glb + call s_mpi_abort('VCFL is greater than 1.0. Exiting ...') + end if + end if + + call s_mpi_barrier() + + end subroutine s_write_run_time_information + + !> The goal of this subroutine is to write to the run-time + !! information file basic footer information applicable to + !! the current computation and to close the file when done. + !! The footer contains the stability criteria extrema over + !! all of the time-steps and the simulation run-time. + subroutine s_close_run_time_information_file + + real(kind(0d0)) :: run_time !< Run-time of the simulation + ! Writing the footer of and closing the run-time information file + write (3, '(A)') '----------------------------------------'// & + '----------------------------------------' + write (3, '(A)') '' + + write (3, '(A,F9.6)') 'ICFL Max: ', icfl_max + if (any(Re_size > 0)) write (3, '(A,F9.6)') 'VCFL Max: ', vcfl_max + if (any(Re_size > 0)) write (3, '(A,F10.6)') 'Rc Min: ', Rc_min + + call cpu_time(run_time) + + write (3, '(A)') '' + write (3, '(A,I0,A)') 'Run-time: ', int(anint(run_time)), 's' + write (3, '(A)') '========================================'// & + '========================================' + close (3) + + end subroutine s_close_run_time_information_file + + subroutine s_finalize_sim_helpers_module() + + ! Deallocating the ICFL, VCFL, CCFL, and Rc stability criteria + @:DEALLOCATE_GLOBAL(icfl_sf) + if (any(Re_size > 0)) then + @:DEALLOCATE_GLOBAL(vcfl_sf, Rc_sf) + end if + + s_read_data_files => null() + s_write_data_files => null() + + end subroutine s_finalize_sim_helpers_module + +end module m_sim_helpers diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 004900eba..05d9d2af0 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -78,15 +78,14 @@ module m_start_up use m_surface_tension use m_body_forces + + use m_sim_helpers ! ========================================================================== implicit none private; public :: s_read_input_file, & s_check_input_file, & - s_read_data_files, & - s_read_serial_data_files, & - s_read_parallel_data_files, & s_populate_grid_variables_buffers, & s_initialize_internal_energy_equations, & s_initialize_modules, s_initialize_gpu_vars, & @@ -94,25 +93,6 @@ module m_start_up s_perform_time_step, s_save_data, & s_save_performance_metrics - abstract interface ! =================================================== - - !! @param q_cons_vf Conservative variables - subroutine s_read_abstract_data_files(q_cons_vf) - - import :: scalar_field, sys_size, pres_field - - type(scalar_field), & - dimension(sys_size), & - intent(inout) :: q_cons_vf - - end subroutine s_read_abstract_data_files - - end interface ! ======================================================== - - type(scalar_field), allocatable, dimension(:) :: grad_x_vf, grad_y_vf, grad_z_vf, norm_vf - - procedure(s_read_abstract_data_files), pointer :: s_read_data_files => null() - contains !> The purpose of this procedure is to first verify that an @@ -232,596 +212,6 @@ contains end subroutine s_check_input_file - !! initial condition and grid data files. The cell-average - !! conservative variables constitute the former, while the - !! cell-boundary locations in x-, y- and z-directions make - !! up the latter. This procedure also calculates the cell- - !! width distributions from the cell-boundary locations. - !! @param q_cons_vf Cell-averaged conservative variables - subroutine s_read_serial_data_files(q_cons_vf) - - type(scalar_field), dimension(sys_size), intent(INOUT) :: q_cons_vf - - character(LEN=path_len + 2*name_len) :: t_step_dir !< - !! Relative path to the starting time-step directory - - character(LEN=path_len + 3*name_len) :: file_path !< - !! Relative path to the grid and conservative variables data files - - logical :: file_exist !< - ! Logical used to check the existence of the data files - - integer :: i, r !< Generic loop iterator - - ! Confirming that the directory from which the initial condition and - ! the grid data files are to be read in exists and exiting otherwise - if (cfl_dt) then - write (t_step_dir, '(A,I0,A,I0)') & - trim(case_dir)//'/p_all/p', proc_rank, '/', n_start - else - write (t_step_dir, '(A,I0,A,I0)') & - trim(case_dir)//'/p_all/p', proc_rank, '/', t_step_start - end if - - file_path = trim(t_step_dir)//'/.' - call my_inquire(file_path, file_exist) - - if (file_exist .neqv. .true.) then - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - - ! Cell-boundary Locations in x-direction =========================== - file_path = trim(t_step_dir)//'/x_cb.dat' - - inquire (FILE=trim(file_path), EXIST=file_exist) - - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) x_cb(-1:m); close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - - dx(0:m) = x_cb(0:m) - x_cb(-1:m - 1) - x_cc(0:m) = x_cb(-1:m - 1) + dx(0:m)/2d0 - - if (ib) then - do i = 1, num_ibs - if (patch_ib(i)%c > 0) then - Np = int((patch_ib(i)%p*patch_ib(i)%c/dx(0))*20) + int(((patch_ib(i)%c - patch_ib(i)%p*patch_ib(i)%c)/dx(0))*20) + 1 - end if - end do - end if - ! ================================================================== - - ! Cell-boundary Locations in y-direction =========================== - if (n > 0) then - - file_path = trim(t_step_dir)//'/y_cb.dat' - - inquire (FILE=trim(file_path), EXIST=file_exist) - - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) y_cb(-1:n); close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - - dy(0:n) = y_cb(0:n) - y_cb(-1:n - 1) - y_cc(0:n) = y_cb(-1:n - 1) + dy(0:n)/2d0 - - end if - ! ================================================================== - - ! Cell-boundary Locations in z-direction =========================== - if (p > 0) then - - file_path = trim(t_step_dir)//'/z_cb.dat' - - inquire (FILE=trim(file_path), EXIST=file_exist) - - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) z_cb(-1:p); close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - - dz(0:p) = z_cb(0:p) - z_cb(-1:p - 1) - z_cc(0:p) = z_cb(-1:p - 1) + dz(0:p)/2d0 - - end if - ! ================================================================== - - do i = 1, sys_size - write (file_path, '(A,I0,A)') & - trim(t_step_dir)//'/q_cons_vf', i, '.dat' - inquire (FILE=trim(file_path), EXIST=file_exist) - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) q_cons_vf(i)%sf(0:m, 0:n, 0:p); close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - end do - - if ((bubbles .eqv. .true.) .or. (hypoelasticity .eqv. .true.)) then - ! Read pb and mv for non-polytropic qbmm - if (qbmm .and. .not. polytropic) then - do i = 1, nb - do r = 1, nnode - write (file_path, '(A,I0,A)') & - trim(t_step_dir)//'/pb', sys_size + (i - 1)*nnode + r, '.dat' - inquire (FILE=trim(file_path), EXIST=file_exist) - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) pb_ts(1)%sf(0:m, 0:n, 0:p, r, i); close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - end do - end do - do i = 1, nb - do r = 1, nnode - write (file_path, '(A,I0,A)') & - trim(t_step_dir)//'/mv', sys_size + (i - 1)*nnode + r, '.dat' - inquire (FILE=trim(file_path), EXIST=file_exist) - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) mv_ts(1)%sf(0:m, 0:n, 0:p, r, i); close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - end do - end do - end if - end if - ! ================================================================== - - ! Read IBM Data ==================================================== - - if (ib) then - do i = 1, num_ibs - write (file_path, '(A,I0,A)') & - trim(t_step_dir)//'/ib.dat' - inquire (FILE=trim(file_path), EXIST=file_exist) - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) ib_markers%sf(0:m, 0:n, 0:p); close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - - if (patch_ib(i)%c > 0) then - - print *, "HERE Np", Np - allocate (airfoil_grid_u(1:Np)) - allocate (airfoil_grid_l(1:Np)) - - write (file_path, '(A)') & - trim(t_step_dir)//'/airfoil_u.dat' - inquire (FILE=trim(file_path), EXIST=file_exist) - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) airfoil_grid_u; close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - - write (file_path, '(A)') & - trim(t_step_dir)//'/airfoil_l.dat' - inquire (FILE=trim(file_path), EXIST=file_exist) - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) airfoil_grid_l; close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - end if - end do - - end if - - end subroutine s_read_serial_data_files - - !! @param q_cons_vf Conservative variables - subroutine s_read_parallel_data_files(q_cons_vf) - - type(scalar_field), & - dimension(sys_size), & - intent(INOUT) :: q_cons_vf - -#ifdef MFC_MPI - - real(kind(0d0)), allocatable, dimension(:) :: x_cb_glb, y_cb_glb, z_cb_glb - - 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) :: 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_start_string - - integer :: i, j - - allocate (x_cb_glb(-1:m_glb)) - allocate (y_cb_glb(-1:n_glb)) - allocate (z_cb_glb(-1:p_glb)) - - ! Read in cell boundary locations in x-direction - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//'x_cb.dat' - inquire (FILE=trim(file_loc), EXIST=file_exist) - - if (file_exist) then - data_size = m_glb + 2 - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - call MPI_FILE_READ(ifile, x_cb_glb, data_size, MPI_DOUBLE_PRECISION, status, ierr) - call MPI_FILE_CLOSE(ifile, ierr) - else - call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') - end if - - ! Assigning local cell boundary locations - x_cb(-1:m) = x_cb_glb((start_idx(1) - 1):(start_idx(1) + m)) - ! Computing the cell width distribution - dx(0:m) = x_cb(0:m) - x_cb(-1:m - 1) - ! Computing the cell center locations - x_cc(0:m) = x_cb(-1:m - 1) + dx(0:m)/2d0 - - if (ib) then - do i = 1, num_ibs - if (patch_ib(i)%c > 0) then - Np = int((patch_ib(i)%p*patch_ib(i)%c/dx(0))*20) + int(((patch_ib(i)%c - patch_ib(i)%p*patch_ib(i)%c)/dx(0))*20) + 1 - allocate (MPI_IO_airfoil_IB_DATA%var(1:2*Np)) - print *, "HERE Np", Np - end if - end do - end if - - if (n > 0) then - ! Read in cell boundary locations in y-direction - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//'y_cb.dat' - inquire (FILE=trim(file_loc), EXIST=file_exist) - - if (file_exist) then - data_size = n_glb + 2 - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - call MPI_FILE_READ(ifile, y_cb_glb, data_size, MPI_DOUBLE_PRECISION, status, ierr) - call MPI_FILE_CLOSE(ifile, ierr) - else - call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') - end if - - ! Assigning local cell boundary locations - y_cb(-1:n) = y_cb_glb((start_idx(2) - 1):(start_idx(2) + n)) - ! Computing the cell width distribution - dy(0:n) = y_cb(0:n) - y_cb(-1:n - 1) - ! Computing the cell center locations - y_cc(0:n) = y_cb(-1:n - 1) + dy(0:n)/2d0 - - if (p > 0) then - ! Read in cell boundary locations in z-direction - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//'z_cb.dat' - inquire (FILE=trim(file_loc), EXIST=file_exist) - - if (file_exist) then - data_size = p_glb + 2 - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - call MPI_FILE_READ(ifile, z_cb_glb, data_size, MPI_DOUBLE_PRECISION, status, ierr) - call MPI_FILE_CLOSE(ifile, ierr) - else - call s_mpi_abort('File '//trim(file_loc)//'is missing. Exiting...') - end if - - ! Assigning local cell boundary locations - z_cb(-1:p) = z_cb_glb((start_idx(3) - 1):(start_idx(3) + p)) - ! Computing the cell width distribution - dz(0:p) = z_cb(0:p) - z_cb(-1:p - 1) - ! Computing the cell center locations - z_cc(0:p) = z_cb(-1:p - 1) + dz(0:p)/2d0 - - end if - end if - - if (file_per_process) then - if (cfl_dt) then - call s_int_to_str(n_start, t_step_start_string) - write (file_loc, '(I0,A1,I7.7,A)') n_start, '_', proc_rank, '.dat' - else - call s_int_to_str(t_step_start, t_step_start_string) - write (file_loc, '(I0,A1,I7.7,A)') t_step_start, '_', proc_rank, '.dat' - end if - - file_loc = trim(case_dir)//'/restart_data/lustre_'//trim(t_step_start_string)//trim(mpiiofs)//trim(file_loc) - inquire (FILE=trim(file_loc), EXIST=file_exist) - - if (file_exist) then - call MPI_FILE_OPEN(MPI_COMM_SELF, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - - ! Initialize MPI data I/O - - if (ib) then - call s_initialize_mpi_data(q_cons_vf, ib_markers) - else - call s_initialize_mpi_data(q_cons_vf) - end if - - ! Size of local arrays - data_size = (m + 1)*(n + 1)*(p + 1) - - ! Resize some integers so MPI can read even the biggest file - m_MOK = int(m_glb + 1, MPI_OFFSET_KIND) - n_MOK = int(n_glb + 1, MPI_OFFSET_KIND) - p_MOK = int(p_glb + 1, MPI_OFFSET_KIND) - WP_MOK = int(8d0, MPI_OFFSET_KIND) - MOK = int(1d0, MPI_OFFSET_KIND) - str_MOK = int(name_len, MPI_OFFSET_KIND) - NVARS_MOK = int(sys_size, MPI_OFFSET_KIND) - - ! Read the data for each variable - if (bubbles .or. hypoelasticity) then - - do i = 1, sys_size!adv_idx%end - var_MOK = int(i, MPI_OFFSET_KIND) - - call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & - MPI_DOUBLE_PRECISION, status, ierr) - end do - !Read pb and mv for non-polytropic qbmm - if (qbmm .and. .not. polytropic) then - do i = sys_size + 1, sys_size + 2*nb*nnode - var_MOK = int(i, MPI_OFFSET_KIND) - - call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & - MPI_DOUBLE_PRECISION, status, ierr) - end do - end if - else - do i = 1, adv_idx%end - var_MOK = int(i, MPI_OFFSET_KIND) - - call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & - MPI_DOUBLE_PRECISION, status, ierr) - end do - end if - - call s_mpi_barrier() - - call MPI_FILE_CLOSE(ifile, ierr) - - if (ib) then - - write (file_loc, '(A)') 'ib.dat' - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) - inquire (FILE=trim(file_loc), EXIST=file_exist) - - if (file_exist) then - - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - - disp = 0 - - call MPI_FILE_SET_VIEW(ifile, disp, MPI_INTEGER, MPI_IO_IB_DATA%view, & - 'native', mpi_info_int, ierr) - call MPI_FILE_READ(ifile, MPI_IO_IB_DATA%var%sf, data_size, & - MPI_INTEGER, status, ierr) - - else - call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') - end if - - end if - - else - call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') - end if - else - ! Open the file to read conservative variables - if (cfl_dt) then - write (file_loc, '(I0,A)') n_start, '.dat' - else - write (file_loc, '(I0,A)') t_step_start, '.dat' - end if - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) - inquire (FILE=trim(file_loc), EXIST=file_exist) - - if (file_exist) then - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - - ! Initialize MPI data I/O - - if (ib) then - call s_initialize_mpi_data(q_cons_vf, ib_markers) - else - - call s_initialize_mpi_data(q_cons_vf) - - end if - - - ! Size of local arrays - data_size = (m + 1)*(n + 1)*(p + 1) - - ! Resize some integers so MPI can read even the biggest file - m_MOK = int(m_glb + 1, MPI_OFFSET_KIND) - n_MOK = int(n_glb + 1, MPI_OFFSET_KIND) - p_MOK = int(p_glb + 1, MPI_OFFSET_KIND) - WP_MOK = int(8d0, MPI_OFFSET_KIND) - MOK = int(1d0, MPI_OFFSET_KIND) - str_MOK = int(name_len, MPI_OFFSET_KIND) - NVARS_MOK = int(sys_size, MPI_OFFSET_KIND) - - ! Read the data for each variable - if (bubbles .or. hypoelasticity) then - - do i = 1, sys_size!adv_idx%end - var_MOK = int(i, MPI_OFFSET_KIND) - ! Initial displacement to skip at beginning of file - disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1) - - call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_DATA%view(i), & - 'native', mpi_info_int, ierr) - call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & - MPI_DOUBLE_PRECISION, status, ierr) - end do - !Read pb and mv for non-polytropic qbmm - if (qbmm .and. .not. polytropic) then - do i = sys_size + 1, sys_size + 2*nb*nnode - var_MOK = int(i, MPI_OFFSET_KIND) - ! Initial displacement to skip at beginning of file - disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1) - - call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_DATA%view(i), & - 'native', mpi_info_int, ierr) - call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & - MPI_DOUBLE_PRECISION, status, ierr) - end do - end if - else - do i = 1, sys_size - var_MOK = int(i, MPI_OFFSET_KIND) - - ! Initial displacement to skip at beginning of file - disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1) - - call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_DATA%view(i), & - 'native', mpi_info_int, ierr) - call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & - MPI_DOUBLE_PRECISION, status, ierr) - - end do - end if - - call s_mpi_barrier() - - call MPI_FILE_CLOSE(ifile, ierr) - - if (ib) then - - write (file_loc, '(A)') 'ib.dat' - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) - inquire (FILE=trim(file_loc), EXIST=file_exist) - - if (file_exist) then - - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - - disp = 0 - - call MPI_FILE_SET_VIEW(ifile, disp, MPI_INTEGER, MPI_IO_IB_DATA%view, & - 'native', mpi_info_int, ierr) - call MPI_FILE_READ(ifile, MPI_IO_IB_DATA%var%sf, data_size, & - MPI_INTEGER, status, ierr) - - else - call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') - end if - - end if - - else - call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') - end if - - end if - - if (ib) then - - do j = 1, num_ibs - if (patch_ib(j)%c > 0) then - - print *, "HERE Np", Np - - allocate (airfoil_grid_u(1:Np)) - allocate (airfoil_grid_l(1:Np)) - - write (file_loc, '(A)') 'airfoil_l.dat' - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) - inquire (FILE=trim(file_loc), EXIST=file_exist) - if (file_exist) then - - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - - ! Initial displacement to skip at beginning of file - disp = 0 - - call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_airfoil_IB_DATA%view(1), & - 'native', mpi_info_int, ierr) - call MPI_FILE_READ(ifile, MPI_IO_airfoil_IB_DATA%var(1:Np), 3*Np, & - MPI_DOUBLE_PRECISION, status, ierr) - - end if - - write (file_loc, '(A)') 'airfoil_u.dat' - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) - inquire (FILE=trim(file_loc), EXIST=file_exist) - if (file_exist) then - - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - - ! Initial displacement to skip at beginning of file - disp = 0 - - call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_airfoil_IB_DATA%view(2), & - 'native', mpi_info_int, ierr) - call MPI_FILE_READ(ifile, MPI_IO_airfoil_IB_DATA%var(Np + 1:2*Np), 3*Np, & - MPI_DOUBLE_PRECISION, status, ierr) - end if - - do i = 1, Np - airfoil_grid_l(i)%x = MPI_IO_airfoil_IB_DATA%var(i)%x - airfoil_grid_l(i)%y = MPI_IO_airfoil_IB_DATA%var(i)%y - end do - - do i = 1, Np - airfoil_grid_u(i)%x = MPI_IO_airfoil_IB_DATA%var(Np + i)%x - airfoil_grid_u(i)%y = MPI_IO_airfoil_IB_DATA%var(Np + i)%y - end do - - end if - end do - end if - - deallocate (x_cb_glb, y_cb_glb, z_cb_glb) - -#endif - - end subroutine s_read_parallel_data_files - !> The purpose of this subroutine is to populate the buffers !! of the grid variables, which are constituted of the cell- !! boundary locations and cell-width distributions, based on @@ -1306,6 +696,8 @@ contains pref = 1d0 end if + call s_initialize_sim_helpers_module() + #if defined(MFC_OpenACC) && defined(MFC_MEMORY_DUMP) call acc_present_dump() #endif @@ -1351,17 +743,8 @@ contains call acc_present_dump() #endif - ! Associate pointers for serial or parallel I/O - if (parallel_io .neqv. .true.) then - s_read_data_files => s_read_serial_data_files - s_write_data_files => s_write_serial_data_files - else - s_read_data_files => s_read_parallel_data_files - s_write_data_files => s_write_parallel_data_files - end if - ! Reading in the user provided initial condition and grid data - call s_read_data_files(q_cons_ts(1)%vf) + call s_read_data_files(q_cons_ts(1)%vf, ib_markers) if (model_eqns == 3) call s_initialize_internal_energy_equations(q_cons_ts(1)%vf) if (ib) call s_ibm_setup() @@ -1497,6 +880,7 @@ contains s_read_data_files => null() s_write_data_files => null() + call s_finalize_sim_helpers_module() call s_finalize_time_steppers_module() call s_finalize_derived_variables_module() call s_finalize_data_output_module() diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index e7327190e..f7a90e6da 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -31,8 +31,6 @@ module m_time_steppers use m_helper - use m_sim_helpers - use m_fftw use m_nvtx @@ -1030,95 +1028,6 @@ contains end subroutine s_apply_bodyforces - subroutine s_comprehensive_debug(q_cons_vf, q_prim_vf, t_step, stage) - - type(scalar_field), dimension(sys_size) :: q_cons_vf, q_prim_vf - integer, intent(in) :: t_step, stage - integer :: j, k, l, i - integer errors - logical :: exists - - character(LEN=name_len) :: file_name = 'comp_debug.txt' - character(LEN=path_len + name_len) :: file_path - character(100) :: str_format - - ! Opening the run-time information file - file_path = trim(case_dir)//'/'//trim(file_name) - - str_format = "(I9, A, I3, A, I4, I4, I4, A, I2, A, I5, A, I5, I5, I5)" - - open (12, FILE=trim(file_path), & - STATUS='replace') - - errors = 0 - - ! Check all variables for NaNs - do i = 1, sys_size - !$acc update host(q_cons_vf(i)%sf) - do l = 0, p - do k = 0, n - do j = 0, m - if (ieee_is_nan(q_cons_vf(i)%sf(j, k, l))) then - write(12, str_format) t_step, " NaN(s) in conservative variables after RK stage ", & - stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & - " (m, n, p)", m, n, p - errors = errors + 1 - end if - end do - end do - end do - end do - - ! Check for invalid volume fractions - do i = advxb, advxe - do l = 0, p - do k = 0, n - do j = 0, m - if (q_cons_vf(i)%sf(j, k, l) < 0d0) then - write(12, str_format) t_step, " Volume fraction < 0 after RK stage ", & - stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & - " (m, n, p)", m, n, p - errors = errors + 1 - elseif (q_cons_vf(i)%sf(j, k, l) > 1d0 + verysmall) then - write(12, str_format) t_step, " Volume fraction > 1 after RK stage ", & - stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & - " (m, n, p)", m, n, p - errors = errors + 1 - end if - end do - end do - end do - end do - - ! Check for invalid densities - do i = contxb, contxe - do l = 0, p - do k = 0, n - do j = 0, m - if (q_cons_vf(advxb + i -1)%sf(j, k, l) < 0d0 .and. q_cons_vf(i)%sf(j, k, l) < 0d0 .or. & - q_cons_vf(advxb + i -1)%sf(j, k, l) > 0d0 .and. q_cons_Vf(i)%sf(j, k, l) < 0d0) then - print*, q_cons_vf(advxb + i - 1)%sf(j, k, l), q_cons_vf(i)%sf(j, k, l) - write(12, str_format) t_step, " Density is negative after RK stage ", & - stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & - " (m, n, p)", m, n, p - errors = errors + 1 - end if - end do - end do - end do - end do - - if (errors /= 0) then - close(12) - call s_write_data_files(q_cons_vf, q_prim_vf, t_step) - call s_mpi_abort("Errors found in conservative variables") - endif - - write(12, "(I3)") -1 - close(12) - - end subroutine s_comprehensive_debug - !> This subroutine saves the temporary q_prim_vf vector !! into the q_prim_ts vector that is then used in p_main !! @param t_step current time-step From dc2a4c426aa977fdc79e97de9552c1bd574da58f Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Wed, 16 Oct 2024 13:40:39 -0400 Subject: [PATCH 06/10] Revert "refactor" This reverts commit f4283a4af8af992fd76c9806bdf1b0ee5883a66f. --- src/simulation/m_data_output.fpp | 291 +++++++ src/simulation/m_sim_helpers.f90 | 250 ++++++ src/simulation/m_sim_helpers.fpp | 1272 ---------------------------- src/simulation/m_start_up.fpp | 628 +++++++++++++- src/simulation/m_time_steppers.fpp | 91 ++ 5 files changed, 1254 insertions(+), 1278 deletions(-) create mode 100644 src/simulation/m_sim_helpers.f90 delete mode 100644 src/simulation/m_sim_helpers.fpp diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 0a48562c1..b5d33da11 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -26,6 +26,8 @@ module m_data_output use m_helper + use m_sim_helpers + use m_delay_file_access use m_ibm @@ -35,11 +37,14 @@ module m_data_output private; public :: s_initialize_data_output_module, & + s_open_run_time_information_file, & s_open_probe_files, & + s_write_run_time_information, & s_write_data_files, & s_write_serial_data_files, & s_write_parallel_data_files, & s_write_probe_files, & + s_close_run_time_information_file, & s_close_probe_files, & s_finalize_data_output_module @@ -65,11 +70,103 @@ module m_data_output end subroutine s_write_abstract_data_files end interface ! ======================================================== +#ifdef CRAY_ACC_WAR + @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), icfl_sf) + @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), vcfl_sf) + @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), ccfl_sf) + @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), Rc_sf) + !$acc declare link(icfl_sf, vcfl_sf, ccfl_sf, Rc_sf) +#else + real(kind(0d0)), allocatable, dimension(:, :, :) :: icfl_sf !< ICFL stability criterion + real(kind(0d0)), allocatable, dimension(:, :, :) :: vcfl_sf !< VCFL stability criterion + real(kind(0d0)), allocatable, dimension(:, :, :) :: ccfl_sf !< CCFL stability criterion + real(kind(0d0)), allocatable, dimension(:, :, :) :: Rc_sf !< Rc stability criterion + !$acc declare create(icfl_sf, vcfl_sf, ccfl_sf, Rc_sf) +#endif + + real(kind(0d0)) :: icfl_max_loc, icfl_max_glb !< ICFL stability extrema on local and global grids + real(kind(0d0)) :: vcfl_max_loc, vcfl_max_glb !< VCFL stability extrema on local and global grids + real(kind(0d0)) :: ccfl_max_loc, ccfl_max_glb !< CCFL stability extrema on local and global grids + real(kind(0d0)) :: Rc_min_loc, Rc_min_glb !< Rc stability extrema on local and global grids + !$acc declare create(icfl_max_loc, icfl_max_glb, vcfl_max_loc, vcfl_max_glb, ccfl_max_loc, ccfl_max_glb, Rc_min_loc, Rc_min_glb) + + !> @name ICFL, VCFL, CCFL and Rc stability criteria extrema over all the time-steps + !> @{ + real(kind(0d0)) :: icfl_max !< ICFL criterion maximum + real(kind(0d0)) :: vcfl_max !< VCFL criterion maximum + real(kind(0d0)) :: ccfl_max !< CCFL criterion maximum + real(kind(0d0)) :: Rc_min !< Rc criterion maximum + !> @} procedure(s_write_abstract_data_files), pointer :: s_write_data_files => null() contains + !> The purpose of this subroutine is to open a new or pre- + !! existing run-time information file and append to it the + !! basic header information relevant to current simulation. + !! In general, this requires generating a table header for + !! those stability criteria which will be written at every + !! time-step. + subroutine s_open_run_time_information_file + + character(LEN=name_len) :: file_name = 'run_time.inf' !< + !! Name of the run-time information file + + character(LEN=path_len + name_len) :: file_path !< + !! Relative path to a file in the case directory + + character(LEN=8) :: file_date !< + !! Creation date of the run-time information file + + ! Opening the run-time information file + file_path = trim(case_dir)//'/'//trim(file_name) + + open (3, FILE=trim(file_path), & + FORM='formatted', & + STATUS='replace') + + write (3, '(A)') 'Description: Stability information at '// & + 'each time-step of the simulation. This' + write (3, '(13X,A)') 'data is composed of the inviscid '// & + 'Courant–Friedrichs–Lewy (ICFL)' + write (3, '(13X,A)') 'number, the viscous CFL (VCFL) number, '// & + 'the capillary CFL (CCFL)' + write (3, '(13X,A)') 'number and the cell Reynolds (Rc) '// & + 'number. Please note that only' + write (3, '(13X,A)') 'those stability conditions pertinent '// & + 'to the physics included in' + write (3, '(13X,A)') 'the current computation are displayed.' + + call date_and_time(DATE=file_date) + + write (3, '(A)') 'Date: '//file_date(5:6)//'/'// & + file_date(7:8)//'/'// & + file_date(3:4) + + write (3, '(A)') ''; write (3, '(A)') '' + + ! Generating table header for the stability criteria to be outputted + if (cfl_dt) then + if (any(Re_size > 0)) then + write (1, '(A)') '==== Time-steps ====== dt ===== Time ======= ICFL '// & + 'Max ==== VCFL Max ====== Rc Min =======' + else + write (1, '(A)') '=========== Time-steps ============== dt ===== Time '// & + '============== ICFL Max =============' + end if + else + if (any(Re_size > 0)) then + write (1, '(A)') '==== Time-steps ====== Time ======= ICFL '// & + 'Max ==== VCFL Max ====== Rc Min =======' + else + write (1, '(A)') '=========== Time-steps ============== Time '// & + '============== ICFL Max =============' + end if + end if + + end subroutine s_open_run_time_information_file + !> This opens a formatted data file where the root processor !! can write out flow probe information subroutine s_open_probe_files @@ -119,6 +216,155 @@ contains end subroutine s_open_probe_files + !> The goal of the procedure is to output to the run-time + !! information file the stability criteria extrema in the + !! entire computational domain and at the given time-step. + !! Moreover, the subroutine is also in charge of tracking + !! these stability criteria extrema over all time-steps. + !! @param q_prim_vf Cell-average primitive variables + !! @param t_step Current time step + subroutine s_write_run_time_information(q_prim_vf, t_step) + + type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf + integer, intent(IN) :: t_step + + real(kind(0d0)), dimension(num_fluids) :: alpha_rho !< Cell-avg. partial density + real(kind(0d0)) :: rho !< Cell-avg. density + real(kind(0d0)), dimension(num_dims) :: vel !< Cell-avg. velocity + real(kind(0d0)) :: vel_sum !< Cell-avg. velocity sum + real(kind(0d0)) :: pres !< Cell-avg. pressure + real(kind(0d0)), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction + real(kind(0d0)) :: gamma !< Cell-avg. sp. heat ratio + real(kind(0d0)) :: pi_inf !< Cell-avg. liquid stiffness function + real(kind(0d0)) :: qv !< Cell-avg. fluid reference energy + real(kind(0d0)) :: c !< Cell-avg. sound speed + real(kind(0d0)) :: E !< Cell-avg. energy + real(kind(0d0)) :: H !< Cell-avg. enthalpy + real(kind(0d0)), dimension(2) :: Re !< Cell-avg. Reynolds numbers + + ! ICFL, VCFL, CCFL and Rc stability criteria extrema for the current + ! time-step and located on both the local (loc) and the global (glb) + ! computational domains + + real(kind(0d0)) :: blkmod1, blkmod2 !< + !! Fluid bulk modulus for Woods mixture sound speed + + integer :: i, j, k, l, q !< Generic loop iterators + + integer :: Nfq + real(kind(0d0)) :: fltr_dtheta !< + !! Modified dtheta accounting for Fourier filtering in azimuthal direction. + + ! Computing Stability Criteria at Current Time-step ================ + !$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho, vel, alpha, Re, fltr_dtheta, Nfq) + do l = 0, p + do k = 0, n + do j = 0, m + + call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l) + + ! Compute mixture sound speed + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, c) + + if (any(Re_size > 0)) then + call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf, vcfl_sf, Rc_sf) + else + call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf) + end if + + end do + end do + end do + ! END: Computing Stability Criteria at Current Time-step =========== + + ! Determining local stability criteria extrema at current time-step + +#ifdef CRAY_ACC_WAR + !$acc update host(icfl_sf) + + if (any(Re_size > 0)) then + !$acc update host(vcfl_sf, Rc_sf) + end if + + icfl_max_loc = maxval(icfl_sf) + + if (any(Re_size > 0)) then + vcfl_max_loc = maxval(vcfl_sf) + Rc_min_loc = minval(Rc_sf) + end if +#else + !$acc kernels + icfl_max_loc = maxval(icfl_sf) + !$acc end kernels + + if (any(Re_size > 0)) then + !$acc kernels + vcfl_max_loc = maxval(vcfl_sf) + Rc_min_loc = minval(Rc_sf) + !$acc end kernels + end if +#endif + + ! Determining global stability criteria extrema at current time-step + if (num_procs > 1) then + call s_mpi_reduce_stability_criteria_extrema(icfl_max_loc, & + vcfl_max_loc, & + ccfl_max_loc, & + Rc_min_loc, & + icfl_max_glb, & + vcfl_max_glb, & + ccfl_max_glb, & + Rc_min_glb) + else + icfl_max_glb = icfl_max_loc + if (any(Re_size > 0)) vcfl_max_glb = vcfl_max_loc + if (any(Re_size > 0)) Rc_min_glb = Rc_min_loc + end if + + ! Determining the stability criteria extrema over all the time-steps + if (icfl_max_glb > icfl_max) icfl_max = icfl_max_glb + + if (any(Re_size > 0)) then + if (vcfl_max_glb > vcfl_max) vcfl_max = vcfl_max_glb + if (Rc_min_glb < Rc_min) Rc_min = Rc_min_glb + end if + + ! Outputting global stability criteria extrema at current time-step + if (proc_rank == 0) then + if (any(Re_size > 0)) then + write (1, '(6X,I8,F10.6,6X,6X,F10.6,6X,F9.6,6X,F9.6,6X,F10.6)') & + t_step, dt, t_step*dt, icfl_max_glb, & + vcfl_max_glb, & + Rc_min_glb + else + write (1, '(13X,I8,14X,F10.6,14X,F10.6,13X,F9.6)') & + t_step, dt, t_step*dt, icfl_max_glb + end if + + if (icfl_max_glb /= icfl_max_glb) then + call s_mpi_abort('ICFL is NaN. Exiting ...') + elseif (icfl_max_glb > 1d0) then + print *, 'icfl', icfl_max_glb + call s_mpi_abort('ICFL is greater than 1.0. Exiting ...') + end if + + do i = chemxb, chemxe + !@:ASSERT(all(q_prim_vf(i)%sf(:,:,:) >= -1d0), "bad conc") + !@:ASSERT(all(q_prim_vf(i)%sf(:,:,:) <= 2d0), "bad conc") + end do + + if (vcfl_max_glb /= vcfl_max_glb) then + call s_mpi_abort('VCFL is NaN. Exiting ...') + elseif (vcfl_max_glb > 1d0) then + print *, 'vcfl', vcfl_max_glb + call s_mpi_abort('VCFL is greater than 1.0. Exiting ...') + end if + end if + + call s_mpi_barrier() + + end subroutine s_write_run_time_information + !> The goal of this subroutine is to output the grid and !! conservative variables data files for given time-step. !! @param q_cons_vf Cell-average conservative variables @@ -1320,6 +1566,33 @@ contains end subroutine s_write_probe_files + !> The goal of this subroutine is to write to the run-time + !! information file basic footer information applicable to + !! the current computation and to close the file when done. + !! The footer contains the stability criteria extrema over + !! all of the time-steps and the simulation run-time. + subroutine s_close_run_time_information_file + + real(kind(0d0)) :: run_time !< Run-time of the simulation + ! Writing the footer of and closing the run-time information file + write (3, '(A)') '----------------------------------------'// & + '----------------------------------------' + write (3, '(A)') '' + + write (3, '(A,F9.6)') 'ICFL Max: ', icfl_max + if (any(Re_size > 0)) write (3, '(A,F9.6)') 'VCFL Max: ', vcfl_max + if (any(Re_size > 0)) write (3, '(A,F10.6)') 'Rc Min: ', Rc_min + + call cpu_time(run_time) + + write (3, '(A)') '' + write (3, '(A,I0,A)') 'Run-time: ', int(anint(run_time)), 's' + write (3, '(A)') '========================================'// & + '========================================' + close (3) + + end subroutine s_close_run_time_information_file + !> Closes probe files subroutine s_close_probe_files @@ -1340,6 +1613,18 @@ contains integer :: i !< Generic loop iterator + ! Allocating/initializing ICFL, VCFL, CCFL and Rc stability criteria + @:ALLOCATE_GLOBAL(icfl_sf(0:m, 0:n, 0:p)) + icfl_max = 0d0 + + if (any(Re_size > 0)) then + @:ALLOCATE_GLOBAL(vcfl_sf(0:m, 0:n, 0:p)) + @:ALLOCATE_GLOBAL(Rc_sf (0:m, 0:n, 0:p)) + + vcfl_max = 0d0 + Rc_min = 1d3 + end if + ! Associating the procedural pointer to the appropriate subroutine ! that will be utilized in the conversion to the mixture variables @@ -1367,6 +1652,12 @@ contains integer :: i !< Generic loop iterator + ! Deallocating the ICFL, VCFL, CCFL, and Rc stability criteria + @:DEALLOCATE_GLOBAL(icfl_sf) + if (any(Re_size > 0)) then + @:DEALLOCATE_GLOBAL(vcfl_sf, Rc_sf) + end if + ! Disassociating the pointer to the procedure that was utilized to ! to convert mixture or species variables to the mixture variables s_convert_to_mixture_variables => null() diff --git a/src/simulation/m_sim_helpers.f90 b/src/simulation/m_sim_helpers.f90 new file mode 100644 index 000000000..be95197f6 --- /dev/null +++ b/src/simulation/m_sim_helpers.f90 @@ -0,0 +1,250 @@ +module m_sim_helpers + + ! Dependencies ============================================================= + use m_derived_types !< Definitions of the derived types + + use m_global_parameters + + use m_variables_conversion + + use m_mpi_proxy + ! ========================================================================== + + implicit none + + private; public :: s_compute_enthalpy, & + s_compute_stability_from_dt, & + s_compute_dt_from_cfl + +contains + + !> Computes enthalpy + !! @param q_prim_vf cell centered primitive variables + !! @param pres mixture pressure + !! @param rho mixture density + !! @param gamma mixture gamma + !! @param pi_inf mixture pi_inf + !! @param Re mixture reynolds number + !! @param H mixture enthalpy + !! @param alpha component alphas + !! @param vel directional velocities + !! @param vel_sum squard sum of velocity components + !! @param j x index + !! @param k y index + !! @param l z index + subroutine s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l) + !$acc routine seq + type(scalar_field), dimension(sys_size) :: q_prim_vf + real(kind(0d0)), dimension(num_fluids) :: alpha_rho + real(kind(0d0)), dimension(num_fluids) :: alpha + real(kind(0d0)), dimension(num_dims) :: vel + real(kind(0d0)) :: rho, gamma, pi_inf, qv, vel_sum, E, H, pres + real(kind(0d0)), dimension(2) :: Re + integer :: i, j, k, l + + do i = 1, num_fluids + alpha_rho(i) = q_prim_vf(i)%sf(j, k, l) + alpha(i) = q_prim_vf(E_idx + i)%sf(j, k, l) + end do + + if (bubbles) then + call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re, j, k, l) + else + call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re, j, k, l) + end if + + do i = 1, num_dims + vel(i) = q_prim_vf(contxe + i)%sf(j, k, l) + end do + + vel_sum = 0d0 + do i = 1, num_dims + vel_sum = vel_sum + vel(i)**2d0 + end do + + pres = q_prim_vf(E_idx)%sf(j, k, l) + + E = gamma*pres + pi_inf + 5d-1*rho*vel_sum + qv + + H = (E + pres)/rho + + end subroutine s_compute_enthalpy + + !> Computes stability criterion for a specified dt + !! @param vel directional velocities + !! @param c mixture speed of sound + !! @param Re_l mixture Reynolds number + !! @param j x index + !! @param k y index + !! @param l z index + !! @param icfl_sf cell centered inviscid cfl number + !! @param vcfl_sf (optional) cell centered viscous cfl number + !! @param Rc_sf (optional) cell centered Rc + subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl_sf, vcfl_sf, Rc_sf) + !$acc routine seq + real(kind(0d0)), dimension(num_dims) :: vel + real(kind(0d0)) :: c, icfl_dt, vcfl_dt, rho + real(kind(0d0)), dimension(0:m, 0:n, 0:p) :: icfl_sf + real(kind(0d0)), dimension(0:m, 0:n, 0:p), optional :: vcfl_sf, Rc_sf + real(kind(0d0)) :: fltr_dtheta !< + !! Modified dtheta accounting for Fourier filtering in azimuthal direction. + integer :: j, k, l + integer :: Nfq + real(kind(0d0)), dimension(2) :: Re_l + + if (grid_geometry == 3) then + if (k == 0) then + fltr_dtheta = 2d0*pi*y_cb(0)/3d0 + elseif (k <= fourier_rings) then + Nfq = min(floor(2d0*real(k, kind(0d0))*pi), (p + 1)/2 + 1) + fltr_dtheta = 2d0*pi*y_cb(k - 1)/real(Nfq, kind(0d0)) + else + fltr_dtheta = y_cb(k - 1)*dz(l) + end if + end if + + if (p > 0) then + !3D + if (grid_geometry == 3) then + icfl_sf(j, k, l) = dt/min(dx(j)/(abs(vel(1)) + c), & + dy(k)/(abs(vel(2)) + c), & + fltr_dtheta/(abs(vel(3)) + c)) + else + icfl_sf(j, k, l) = dt/min(dx(j)/(abs(vel(1)) + c), & + dy(k)/(abs(vel(2)) + c), & + dz(l)/(abs(vel(3)) + c)) + end if + + if (any(Re_size > 0)) then + + if (grid_geometry == 3) then + vcfl_sf(j, k, l) = maxval(dt/Re_l/rho) & + /min(dx(j), dy(k), fltr_dtheta)**2d0 + + Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), & + dy(k)*(abs(vel(2)) + c), & + fltr_dtheta*(abs(vel(3)) + c)) & + /maxval(1d0/Re_l) + else + vcfl_sf(j, k, l) = maxval(dt/Re_l/rho) & + /min(dx(j), dy(k), dz(l))**2d0 + + Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), & + dy(k)*(abs(vel(2)) + c), & + dz(l)*(abs(vel(3)) + c)) & + /maxval(1d0/Re_l) + end if + + end if + + elseif (n > 0) then + !2D + icfl_sf(j, k, l) = dt/min(dx(j)/(abs(vel(1)) + c), & + dy(k)/(abs(vel(2)) + c)) + + if (any(Re_size > 0)) then + + vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/min(dx(j), dy(k))**2d0 + + Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), & + dy(k)*(abs(vel(2)) + c)) & + /maxval(1d0/Re_l) + + end if + + else + !1D + icfl_sf(j, k, l) = (dt/dx(j))*(abs(vel(1)) + c) + + if (any(Re_size > 0)) then + + vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/dx(j)**2d0 + + Rc_sf(j, k, l) = dx(j)*(abs(vel(1)) + c)/maxval(1d0/Re_l) + + end if + + end if + + end subroutine s_compute_stability_from_dt + + !> Computes dt for a specified CFL number + !! @param vel directional velocities + !! @param max_dt cell centered maximum dt + !! @param rho cell centered density + !! @param Re_l cell centered Reynolds number + !! @param j x coordinate + !! @param k y coordinate + !! @param l z coordinate + subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l) + !$acc routine seq + real(kind(0d0)), dimension(num_dims) :: vel + real(kind(0d0)) :: c, icfl_dt, vcfl_dt, rho + real(kind(0d0)), dimension(0:m, 0:n, 0:p) :: max_dt + real(kind(0d0)) :: fltr_dtheta !< + !! Modified dtheta accounting for Fourier filtering in azimuthal direction. + integer :: j, k, l + integer :: Nfq + real(kind(0d0)), dimension(2) :: Re_l + + if (grid_geometry == 3) then + if (k == 0) then + fltr_dtheta = 2d0*pi*y_cb(0)/3d0 + elseif (k <= fourier_rings) then + Nfq = min(floor(2d0*real(k, kind(0d0))*pi), (p + 1)/2 + 1) + fltr_dtheta = 2d0*pi*y_cb(k - 1)/real(Nfq, kind(0d0)) + else + fltr_dtheta = y_cb(k - 1)*dz(l) + end if + end if + + if (p > 0) then + !3D + if (grid_geometry == 3) then + icfl_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), & + dy(k)/(abs(vel(2)) + c), & + fltr_dtheta/(abs(vel(3)) + c)) + else + icfl_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), & + dy(k)/(abs(vel(2)) + c), & + dz(l)/(abs(vel(3)) + c)) + end if + + if (any(Re_size > 0)) then + if (grid_geometry == 3) then + vcfl_dt = cfl_target*(min(dx(j), dy(k), fltr_dtheta)**2d0) & + /minval(1/(rho*Re_l)) + else + vcfl_dt = cfl_target*(min(dx(j), dy(k), dz(l))**2d0) & + /minval(1/(rho*Re_l)) + end if + end if + + elseif (n > 0) then + !2D + icfl_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), & + dy(k)/(abs(vel(2)) + c)) + + if (any(Re_size > 0)) then + vcfl_dt = cfl_target*(min(dx(j), dy(k))**2d0)/maxval((1/Re_l)/rho) + end if + + else + !1D + icfl_dt = cfl_target*(dx(j)/(abs(vel(1)) + c)) + + if (any(Re_size > 0)) then + vcfl_dt = cfl_target*(dx(j)**2d0)/minval(1/(rho*Re_l)) + end if + + end if + + if (any(re_size > 0)) then + max_dt(j, k, l) = min(icfl_dt, vcfl_dt) + else + max_dt(j, k, l) = icfl_dt + end if + + end subroutine s_compute_dt_from_cfl + +end module m_sim_helpers diff --git a/src/simulation/m_sim_helpers.fpp b/src/simulation/m_sim_helpers.fpp deleted file mode 100644 index f533ff266..000000000 --- a/src/simulation/m_sim_helpers.fpp +++ /dev/null @@ -1,1272 +0,0 @@ -#:include 'macros.fpp' - -module m_sim_helpers - - ! Dependencies ============================================================= - use m_derived_types !< Definitions of the derived types - - use m_global_parameters - - use m_variables_conversion - - use m_mpi_proxy - - use m_data_output - - use m_compile_specific - ! ========================================================================== - - implicit none - - private; public :: s_compute_enthalpy, & - s_read_data_files, & - s_read_serial_data_files, & - s_read_parallel_data_files, & - s_compute_stability_from_dt, & - s_compute_dt_from_cfl, & - s_comprehensive_debug, & - s_open_run_time_information_file, & - s_write_run_time_information, & - s_close_run_time_information_file, & - s_initialize_sim_helpers_module, & - s_finalize_sim_helpers_module - - abstract interface ! =================================================== - - !! @param q_cons_vf Conservative variables - subroutine s_read_abstract_data_files(q_cons_vf, ib_markers) - - import :: scalar_field, integer_field, sys_size, pres_field - - type(scalar_field), & - dimension(sys_size), & - intent(inout) :: q_cons_vf - - type(integer_field) :: ib_markers - - end subroutine s_read_abstract_data_files - - end interface ! ======================================================== - - procedure(s_read_abstract_data_files), pointer :: s_read_data_files => null() - -#ifdef CRAY_ACC_WAR - @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), icfl_sf) - @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), vcfl_sf) - @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), ccfl_sf) - @:CRAY_DECLARE_GLOBAL(real(kind(0d0)), dimension(:, :, :), Rc_sf) - !$acc declare link(icfl_sf, vcfl_sf, ccfl_sf, Rc_sf) -#else - real(kind(0d0)), allocatable, dimension(:, :, :) :: icfl_sf !< ICFL stability criterion - real(kind(0d0)), allocatable, dimension(:, :, :) :: vcfl_sf !< VCFL stability criterion - real(kind(0d0)), allocatable, dimension(:, :, :) :: ccfl_sf !< CCFL stability criterion - real(kind(0d0)), allocatable, dimension(:, :, :) :: Rc_sf !< Rc stability criterion - !$acc declare create(icfl_sf, vcfl_sf, ccfl_sf, Rc_sf) -#endif - - real(kind(0d0)) :: icfl_max_loc, icfl_max_glb !< ICFL stability extrema on local and global grids - real(kind(0d0)) :: vcfl_max_loc, vcfl_max_glb !< VCFL stability extrema on local and global grids - real(kind(0d0)) :: ccfl_max_loc, ccfl_max_glb !< CCFL stability extrema on local and global grids - real(kind(0d0)) :: Rc_min_loc, Rc_min_glb !< Rc stability extrema on local and global grids - !$acc declare create(icfl_max_loc, icfl_max_glb, vcfl_max_loc, vcfl_max_glb, ccfl_max_loc, ccfl_max_glb, Rc_min_loc, Rc_min_glb) - - !> @name ICFL, VCFL, CCFL and Rc stability criteria extrema over all the time-steps - !> @{ - real(kind(0d0)) :: icfl_max !< ICFL criterion maximum - real(kind(0d0)) :: vcfl_max !< VCFL criterion maximum - real(kind(0d0)) :: ccfl_max !< CCFL criterion maximum - real(kind(0d0)) :: Rc_min !< Rc criterion maximum - !> @} - -contains - - subroutine s_initialize_sim_helpers_module() - - ! Allocating/initializing ICFL, VCFL, CCFL and Rc stability criteria - @:ALLOCATE_GLOBAL(icfl_sf(0:m, 0:n, 0:p)) - icfl_max = 0d0 - - if (any(Re_size > 0)) then - @:ALLOCATE_GLOBAL(vcfl_sf(0:m, 0:n, 0:p)) - @:ALLOCATE_GLOBAL(Rc_sf (0:m, 0:n, 0:p)) - - vcfl_max = 0d0 - Rc_min = 1d3 - end if - - ! Associate pointers for serial or parallel I/O - if (parallel_io .neqv. .true.) then - s_read_data_files => s_read_serial_data_files - s_write_data_files => s_write_serial_data_files - else - s_read_data_files => s_read_parallel_data_files - s_write_data_files => s_write_parallel_data_files - end if - - end subroutine s_initialize_sim_helpers_module - - !> Computes enthalpy - !! @param q_prim_vf cell centered primitive variables - !! @param pres mixture pressure - !! @param rho mixture density - !! @param gamma mixture gamma - !! @param pi_inf mixture pi_inf - !! @param Re mixture reynolds number - !! @param H mixture enthalpy - !! @param alpha component alphas - !! @param vel directional velocities - !! @param vel_sum squard sum of velocity components - !! @param j x index - !! @param k y index - !! @param l z index - subroutine s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l) - !$acc routine seq - type(scalar_field), dimension(sys_size) :: q_prim_vf - real(kind(0d0)), dimension(num_fluids) :: alpha_rho - real(kind(0d0)), dimension(num_fluids) :: alpha - real(kind(0d0)), dimension(num_dims) :: vel - real(kind(0d0)) :: rho, gamma, pi_inf, qv, vel_sum, E, H, pres - real(kind(0d0)), dimension(2) :: Re - integer :: i, j, k, l - - do i = 1, num_fluids - alpha_rho(i) = q_prim_vf(i)%sf(j, k, l) - alpha(i) = q_prim_vf(E_idx + i)%sf(j, k, l) - end do - - if (bubbles) then - call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re, j, k, l) - else - call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re, j, k, l) - end if - - do i = 1, num_dims - vel(i) = q_prim_vf(contxe + i)%sf(j, k, l) - end do - - vel_sum = 0d0 - do i = 1, num_dims - vel_sum = vel_sum + vel(i)**2d0 - end do - - pres = q_prim_vf(E_idx)%sf(j, k, l) - - E = gamma*pres + pi_inf + 5d-1*rho*vel_sum + qv - - H = (E + pres)/rho - - end subroutine s_compute_enthalpy - - !> Computes stability criterion for a specified dt - !! @param vel directional velocities - !! @param c mixture speed of sound - !! @param Re_l mixture Reynolds number - !! @param j x index - !! @param k y index - !! @param l z index - !! @param icfl_sf cell centered inviscid cfl number - !! @param vcfl_sf (optional) cell centered viscous cfl number - !! @param Rc_sf (optional) cell centered Rc - subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl_sf, vcfl_sf, Rc_sf) - !$acc routine seq - real(kind(0d0)), dimension(num_dims) :: vel - real(kind(0d0)) :: c, icfl_dt, vcfl_dt, rho - real(kind(0d0)), dimension(0:m, 0:n, 0:p) :: icfl_sf - real(kind(0d0)), dimension(0:m, 0:n, 0:p), optional :: vcfl_sf, Rc_sf - real(kind(0d0)) :: fltr_dtheta !< - !! Modified dtheta accounting for Fourier filtering in azimuthal direction. - integer :: j, k, l - integer :: Nfq - real(kind(0d0)), dimension(2) :: Re_l - - if (grid_geometry == 3) then - if (k == 0) then - fltr_dtheta = 2d0*pi*y_cb(0)/3d0 - elseif (k <= fourier_rings) then - Nfq = min(floor(2d0*real(k, kind(0d0))*pi), (p + 1)/2 + 1) - fltr_dtheta = 2d0*pi*y_cb(k - 1)/real(Nfq, kind(0d0)) - else - fltr_dtheta = y_cb(k - 1)*dz(l) - end if - end if - - if (p > 0) then - !3D - if (grid_geometry == 3) then - icfl_sf(j, k, l) = dt/min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c), & - fltr_dtheta/(abs(vel(3)) + c)) - else - icfl_sf(j, k, l) = dt/min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c), & - dz(l)/(abs(vel(3)) + c)) - end if - - if (any(Re_size > 0)) then - - if (grid_geometry == 3) then - vcfl_sf(j, k, l) = maxval(dt/Re_l/rho) & - /min(dx(j), dy(k), fltr_dtheta)**2d0 - - Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), & - dy(k)*(abs(vel(2)) + c), & - fltr_dtheta*(abs(vel(3)) + c)) & - /maxval(1d0/Re_l) - else - vcfl_sf(j, k, l) = maxval(dt/Re_l/rho) & - /min(dx(j), dy(k), dz(l))**2d0 - - Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), & - dy(k)*(abs(vel(2)) + c), & - dz(l)*(abs(vel(3)) + c)) & - /maxval(1d0/Re_l) - end if - - end if - - elseif (n > 0) then - !2D - icfl_sf(j, k, l) = dt/min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c)) - - if (any(Re_size > 0)) then - - vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/min(dx(j), dy(k))**2d0 - - Rc_sf(j, k, l) = min(dx(j)*(abs(vel(1)) + c), & - dy(k)*(abs(vel(2)) + c)) & - /maxval(1d0/Re_l) - - end if - - else - !1D - icfl_sf(j, k, l) = (dt/dx(j))*(abs(vel(1)) + c) - - if (any(Re_size > 0)) then - - vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/dx(j)**2d0 - - Rc_sf(j, k, l) = dx(j)*(abs(vel(1)) + c)/maxval(1d0/Re_l) - - end if - - end if - - end subroutine s_compute_stability_from_dt - - !> Computes dt for a specified CFL number - !! @param vel directional velocities - !! @param max_dt cell centered maximum dt - !! @param rho cell centered density - !! @param Re_l cell centered Reynolds number - !! @param j x coordinate - !! @param k y coordinate - !! @param l z coordinate - subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l) - !$acc routine seq - real(kind(0d0)), dimension(num_dims) :: vel - real(kind(0d0)) :: c, icfl_dt, vcfl_dt, rho - real(kind(0d0)), dimension(0:m, 0:n, 0:p) :: max_dt - real(kind(0d0)) :: fltr_dtheta !< - !! Modified dtheta accounting for Fourier filtering in azimuthal direction. - integer :: j, k, l - integer :: Nfq - real(kind(0d0)), dimension(2) :: Re_l - - if (grid_geometry == 3) then - if (k == 0) then - fltr_dtheta = 2d0*pi*y_cb(0)/3d0 - elseif (k <= fourier_rings) then - Nfq = min(floor(2d0*real(k, kind(0d0))*pi), (p + 1)/2 + 1) - fltr_dtheta = 2d0*pi*y_cb(k - 1)/real(Nfq, kind(0d0)) - else - fltr_dtheta = y_cb(k - 1)*dz(l) - end if - end if - - if (p > 0) then - !3D - if (grid_geometry == 3) then - icfl_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c), & - fltr_dtheta/(abs(vel(3)) + c)) - else - icfl_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c), & - dz(l)/(abs(vel(3)) + c)) - end if - - if (any(Re_size > 0)) then - if (grid_geometry == 3) then - vcfl_dt = cfl_target*(min(dx(j), dy(k), fltr_dtheta)**2d0) & - /minval(1/(rho*Re_l)) - else - vcfl_dt = cfl_target*(min(dx(j), dy(k), dz(l))**2d0) & - /minval(1/(rho*Re_l)) - end if - end if - - elseif (n > 0) then - !2D - icfl_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), & - dy(k)/(abs(vel(2)) + c)) - - if (any(Re_size > 0)) then - vcfl_dt = cfl_target*(min(dx(j), dy(k))**2d0)/maxval((1/Re_l)/rho) - end if - - else - !1D - icfl_dt = cfl_target*(dx(j)/(abs(vel(1)) + c)) - - if (any(Re_size > 0)) then - vcfl_dt = cfl_target*(dx(j)**2d0)/minval(1/(rho*Re_l)) - end if - - end if - - if (any(re_size > 0)) then - max_dt(j, k, l) = min(icfl_dt, vcfl_dt) - else - max_dt(j, k, l) = icfl_dt - end if - - end subroutine s_compute_dt_from_cfl - - subroutine s_comprehensive_debug(q_cons_vf, q_prim_vf, t_step, stage) - - type(scalar_field), dimension(sys_size) :: q_cons_vf, q_prim_vf - integer, intent(in) :: t_step, stage - integer :: j, k, l, i - integer errors - logical :: exists - - character(LEN=name_len) :: file_name = 'comp_debug.txt' - character(LEN=path_len + name_len) :: file_path - character(100) :: str_format - - ! Opening the run-time information file - file_path = trim(case_dir)//'/'//trim(file_name) - - str_format = "(I9, A, I3, A, I4, I4, I4, A, I2, A, I5, A, I5, I5, I5)" - - open (12, FILE=trim(file_path), & - STATUS='replace') - - errors = 0 - - ! Check all variables for NaNs - do i = 1, sys_size - !$acc update host(q_cons_vf(i)%sf) - do l = 0, p - do k = 0, n - do j = 0, m - if (ieee_is_nan(q_cons_vf(i)%sf(j, k, l))) then - write(12, str_format) t_step, " NaN(s) in conservative variables after RK stage ", & - stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & - " (m, n, p)", m, n, p - errors = errors + 1 - end if - end do - end do - end do - end do - - ! Check for invalid volume fractions - do i = advxb, advxe - do l = 0, p - do k = 0, n - do j = 0, m - if (q_cons_vf(i)%sf(j, k, l) < 0d0) then - write(12, str_format) t_step, " Volume fraction < 0 after RK stage ", & - stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & - " (m, n, p)", m, n, p - errors = errors + 1 - elseif (q_cons_vf(i)%sf(j, k, l) > 1d0 + verysmall) then - write(12, str_format) t_step, " Volume fraction > 1 after RK stage ", & - stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & - " (m, n, p)", m, n, p - errors = errors + 1 - end if - end do - end do - end do - end do - - ! Check for invalid densities - do i = contxb, contxe - do l = 0, p - do k = 0, n - do j = 0, m - if (q_cons_vf(advxb + i -1)%sf(j, k, l) < 0d0 .and. q_cons_vf(i)%sf(j, k, l) < 0d0 .or. & - q_cons_vf(advxb + i -1)%sf(j, k, l) > 0d0 .and. q_cons_Vf(i)%sf(j, k, l) < 0d0) then - print*, q_cons_vf(advxb + i - 1)%sf(j, k, l), q_cons_vf(i)%sf(j, k, l) - write(12, str_format) t_step, " Density is negative after RK stage ", & - stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & - " (m, n, p)", m, n, p - errors = errors + 1 - end if - end do - end do - end do - end do - - if (errors /= 0) then - close(12) - call s_write_data_files(q_cons_vf, q_prim_vf, t_step) - call s_mpi_abort("Errors found in conservative variables") - endif - - write(12, "(I3)") -1 - close(12) - - end subroutine s_comprehensive_debug - - !! initial condition and grid data files. The cell-average - !! conservative variables constitute the former, while the - !! cell-boundary locations in x-, y- and z-directions make - !! up the latter. This procedure also calculates the cell- - !! width distributions from the cell-boundary locations. - !! @param q_cons_vf Cell-averaged conservative variables - subroutine s_read_serial_data_files(q_cons_vf, ib_markers) - - type(scalar_field), dimension(sys_size), intent(INOUT) :: q_cons_vf - type(integer_field) :: ib_markers - - character(LEN=path_len + 2*name_len) :: t_step_dir !< - !! Relative path to the starting time-step directory - - character(LEN=path_len + 3*name_len) :: file_path !< - !! Relative path to the grid and conservative variables data files - - logical :: file_exist !< - ! Logical used to check the existence of the data files - - integer :: i, r !< Generic loop iterator - - ! Confirming that the directory from which the initial condition and - ! the grid data files are to be read in exists and exiting otherwise - if (cfl_dt) then - write (t_step_dir, '(A,I0,A,I0)') & - trim(case_dir)//'/p_all/p', proc_rank, '/', n_start - else - write (t_step_dir, '(A,I0,A,I0)') & - trim(case_dir)//'/p_all/p', proc_rank, '/', t_step_start - end if - - file_path = trim(t_step_dir)//'/.' - call my_inquire(file_path, file_exist) - - if (file_exist .neqv. .true.) then - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - - ! Cell-boundary Locations in x-direction =========================== - file_path = trim(t_step_dir)//'/x_cb.dat' - - inquire (FILE=trim(file_path), EXIST=file_exist) - - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) x_cb(-1:m); close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - - dx(0:m) = x_cb(0:m) - x_cb(-1:m - 1) - x_cc(0:m) = x_cb(-1:m - 1) + dx(0:m)/2d0 - - if (ib) then - do i = 1, num_ibs - if (patch_ib(i)%c > 0) then - Np = int((patch_ib(i)%p*patch_ib(i)%c/dx(0))*20) + int(((patch_ib(i)%c - patch_ib(i)%p*patch_ib(i)%c)/dx(0))*20) + 1 - end if - end do - end if - ! ================================================================== - - ! Cell-boundary Locations in y-direction =========================== - if (n > 0) then - - file_path = trim(t_step_dir)//'/y_cb.dat' - - inquire (FILE=trim(file_path), EXIST=file_exist) - - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) y_cb(-1:n); close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - - dy(0:n) = y_cb(0:n) - y_cb(-1:n - 1) - y_cc(0:n) = y_cb(-1:n - 1) + dy(0:n)/2d0 - - end if - ! ================================================================== - - ! Cell-boundary Locations in z-direction =========================== - if (p > 0) then - - file_path = trim(t_step_dir)//'/z_cb.dat' - - inquire (FILE=trim(file_path), EXIST=file_exist) - - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) z_cb(-1:p); close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - - dz(0:p) = z_cb(0:p) - z_cb(-1:p - 1) - z_cc(0:p) = z_cb(-1:p - 1) + dz(0:p)/2d0 - - end if - ! ================================================================== - - do i = 1, sys_size - write (file_path, '(A,I0,A)') & - trim(t_step_dir)//'/q_cons_vf', i, '.dat' - inquire (FILE=trim(file_path), EXIST=file_exist) - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) q_cons_vf(i)%sf(0:m, 0:n, 0:p); close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - end do - - if ((bubbles .eqv. .true.) .or. (hypoelasticity .eqv. .true.)) then - ! Read pb and mv for non-polytropic qbmm - if (qbmm .and. .not. polytropic) then - do i = 1, nb - do r = 1, nnode - write (file_path, '(A,I0,A)') & - trim(t_step_dir)//'/pb', sys_size + (i - 1)*nnode + r, '.dat' - inquire (FILE=trim(file_path), EXIST=file_exist) - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) pb_ts(1)%sf(0:m, 0:n, 0:p, r, i); close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - end do - end do - do i = 1, nb - do r = 1, nnode - write (file_path, '(A,I0,A)') & - trim(t_step_dir)//'/mv', sys_size + (i - 1)*nnode + r, '.dat' - inquire (FILE=trim(file_path), EXIST=file_exist) - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) mv_ts(1)%sf(0:m, 0:n, 0:p, r, i); close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - end do - end do - end if - end if - ! ================================================================== - - ! Read IBM Data ==================================================== - - if (ib) then - do i = 1, num_ibs - write (file_path, '(A,I0,A)') & - trim(t_step_dir)//'/ib.dat' - inquire (FILE=trim(file_path), EXIST=file_exist) - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) ib_markers%sf(0:m, 0:n, 0:p); close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - - if (patch_ib(i)%c > 0) then - - print *, "HERE Np", Np - allocate (airfoil_grid_u(1:Np)) - allocate (airfoil_grid_l(1:Np)) - - write (file_path, '(A)') & - trim(t_step_dir)//'/airfoil_u.dat' - inquire (FILE=trim(file_path), EXIST=file_exist) - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) airfoil_grid_u; close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - - write (file_path, '(A)') & - trim(t_step_dir)//'/airfoil_l.dat' - inquire (FILE=trim(file_path), EXIST=file_exist) - if (file_exist) then - open (2, FILE=trim(file_path), & - FORM='unformatted', & - ACTION='read', & - STATUS='old') - read (2) airfoil_grid_l; close (2) - else - call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') - end if - end if - end do - - end if - - end subroutine s_read_serial_data_files - - !! @param q_cons_vf Conservative variables - subroutine s_read_parallel_data_files(q_cons_vf, ib_markers) - - type(scalar_field), & - dimension(sys_size), & - intent(INOUT) :: q_cons_vf - type(integer_field) :: ib_markers - -#ifdef MFC_MPI - - real(kind(0d0)), allocatable, dimension(:) :: x_cb_glb, y_cb_glb, z_cb_glb - - 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) :: 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_start_string - - integer :: i, j - - allocate (x_cb_glb(-1:m_glb)) - allocate (y_cb_glb(-1:n_glb)) - allocate (z_cb_glb(-1:p_glb)) - - ! Read in cell boundary locations in x-direction - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//'x_cb.dat' - inquire (FILE=trim(file_loc), EXIST=file_exist) - - if (file_exist) then - data_size = m_glb + 2 - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - call MPI_FILE_READ(ifile, x_cb_glb, data_size, MPI_DOUBLE_PRECISION, status, ierr) - call MPI_FILE_CLOSE(ifile, ierr) - else - call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') - end if - - ! Assigning local cell boundary locations - x_cb(-1:m) = x_cb_glb((start_idx(1) - 1):(start_idx(1) + m)) - ! Computing the cell width distribution - dx(0:m) = x_cb(0:m) - x_cb(-1:m - 1) - ! Computing the cell center locations - x_cc(0:m) = x_cb(-1:m - 1) + dx(0:m)/2d0 - - if (ib) then - do i = 1, num_ibs - if (patch_ib(i)%c > 0) then - Np = int((patch_ib(i)%p*patch_ib(i)%c/dx(0))*20) + int(((patch_ib(i)%c - patch_ib(i)%p*patch_ib(i)%c)/dx(0))*20) + 1 - allocate (MPI_IO_airfoil_IB_DATA%var(1:2*Np)) - print *, "HERE Np", Np - end if - end do - end if - - if (n > 0) then - ! Read in cell boundary locations in y-direction - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//'y_cb.dat' - inquire (FILE=trim(file_loc), EXIST=file_exist) - - if (file_exist) then - data_size = n_glb + 2 - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - call MPI_FILE_READ(ifile, y_cb_glb, data_size, MPI_DOUBLE_PRECISION, status, ierr) - call MPI_FILE_CLOSE(ifile, ierr) - else - call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') - end if - - ! Assigning local cell boundary locations - y_cb(-1:n) = y_cb_glb((start_idx(2) - 1):(start_idx(2) + n)) - ! Computing the cell width distribution - dy(0:n) = y_cb(0:n) - y_cb(-1:n - 1) - ! Computing the cell center locations - y_cc(0:n) = y_cb(-1:n - 1) + dy(0:n)/2d0 - - if (p > 0) then - ! Read in cell boundary locations in z-direction - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//'z_cb.dat' - inquire (FILE=trim(file_loc), EXIST=file_exist) - - if (file_exist) then - data_size = p_glb + 2 - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - call MPI_FILE_READ(ifile, z_cb_glb, data_size, MPI_DOUBLE_PRECISION, status, ierr) - call MPI_FILE_CLOSE(ifile, ierr) - else - call s_mpi_abort('File '//trim(file_loc)//'is missing. Exiting...') - end if - - ! Assigning local cell boundary locations - z_cb(-1:p) = z_cb_glb((start_idx(3) - 1):(start_idx(3) + p)) - ! Computing the cell width distribution - dz(0:p) = z_cb(0:p) - z_cb(-1:p - 1) - ! Computing the cell center locations - z_cc(0:p) = z_cb(-1:p - 1) + dz(0:p)/2d0 - - end if - end if - - if (file_per_process) then - if (cfl_dt) then - call s_int_to_str(n_start, t_step_start_string) - write (file_loc, '(I0,A1,I7.7,A)') n_start, '_', proc_rank, '.dat' - else - call s_int_to_str(t_step_start, t_step_start_string) - write (file_loc, '(I0,A1,I7.7,A)') t_step_start, '_', proc_rank, '.dat' - end if - - file_loc = trim(case_dir)//'/restart_data/lustre_'//trim(t_step_start_string)//trim(mpiiofs)//trim(file_loc) - inquire (FILE=trim(file_loc), EXIST=file_exist) - - if (file_exist) then - call MPI_FILE_OPEN(MPI_COMM_SELF, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - - ! Initialize MPI data I/O - - if (ib) then - call s_initialize_mpi_data(q_cons_vf, ib_markers) - else - call s_initialize_mpi_data(q_cons_vf) - end if - - ! Size of local arrays - data_size = (m + 1)*(n + 1)*(p + 1) - - ! Resize some integers so MPI can read even the biggest file - m_MOK = int(m_glb + 1, MPI_OFFSET_KIND) - n_MOK = int(n_glb + 1, MPI_OFFSET_KIND) - p_MOK = int(p_glb + 1, MPI_OFFSET_KIND) - WP_MOK = int(8d0, MPI_OFFSET_KIND) - MOK = int(1d0, MPI_OFFSET_KIND) - str_MOK = int(name_len, MPI_OFFSET_KIND) - NVARS_MOK = int(sys_size, MPI_OFFSET_KIND) - - ! Read the data for each variable - if (bubbles .or. hypoelasticity) then - - do i = 1, sys_size!adv_idx%end - var_MOK = int(i, MPI_OFFSET_KIND) - - call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & - MPI_DOUBLE_PRECISION, status, ierr) - end do - !Read pb and mv for non-polytropic qbmm - if (qbmm .and. .not. polytropic) then - do i = sys_size + 1, sys_size + 2*nb*nnode - var_MOK = int(i, MPI_OFFSET_KIND) - - call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & - MPI_DOUBLE_PRECISION, status, ierr) - end do - end if - else - do i = 1, adv_idx%end - var_MOK = int(i, MPI_OFFSET_KIND) - - call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & - MPI_DOUBLE_PRECISION, status, ierr) - end do - end if - - call s_mpi_barrier() - - call MPI_FILE_CLOSE(ifile, ierr) - - if (ib) then - - write (file_loc, '(A)') 'ib.dat' - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) - inquire (FILE=trim(file_loc), EXIST=file_exist) - - if (file_exist) then - - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - - disp = 0 - - call MPI_FILE_SET_VIEW(ifile, disp, MPI_INTEGER, MPI_IO_IB_DATA%view, & - 'native', mpi_info_int, ierr) - call MPI_FILE_READ(ifile, MPI_IO_IB_DATA%var%sf, data_size, & - MPI_INTEGER, status, ierr) - - else - call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') - end if - - end if - - else - call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') - end if - else - ! Open the file to read conservative variables - if (cfl_dt) then - write (file_loc, '(I0,A)') n_start, '.dat' - else - write (file_loc, '(I0,A)') t_step_start, '.dat' - end if - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) - inquire (FILE=trim(file_loc), EXIST=file_exist) - - if (file_exist) then - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - - ! Initialize MPI data I/O - - if (ib) then - call s_initialize_mpi_data(q_cons_vf, ib_markers) - else - - call s_initialize_mpi_data(q_cons_vf) - - end if - - - ! Size of local arrays - data_size = (m + 1)*(n + 1)*(p + 1) - - ! Resize some integers so MPI can read even the biggest file - m_MOK = int(m_glb + 1, MPI_OFFSET_KIND) - n_MOK = int(n_glb + 1, MPI_OFFSET_KIND) - p_MOK = int(p_glb + 1, MPI_OFFSET_KIND) - WP_MOK = int(8d0, MPI_OFFSET_KIND) - MOK = int(1d0, MPI_OFFSET_KIND) - str_MOK = int(name_len, MPI_OFFSET_KIND) - NVARS_MOK = int(sys_size, MPI_OFFSET_KIND) - - ! Read the data for each variable - if (bubbles .or. hypoelasticity) then - - do i = 1, sys_size!adv_idx%end - var_MOK = int(i, MPI_OFFSET_KIND) - ! Initial displacement to skip at beginning of file - disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1) - - call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_DATA%view(i), & - 'native', mpi_info_int, ierr) - call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & - MPI_DOUBLE_PRECISION, status, ierr) - end do - !Read pb and mv for non-polytropic qbmm - if (qbmm .and. .not. polytropic) then - do i = sys_size + 1, sys_size + 2*nb*nnode - var_MOK = int(i, MPI_OFFSET_KIND) - ! Initial displacement to skip at beginning of file - disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1) - - call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_DATA%view(i), & - 'native', mpi_info_int, ierr) - call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & - MPI_DOUBLE_PRECISION, status, ierr) - end do - end if - else - do i = 1, sys_size - var_MOK = int(i, MPI_OFFSET_KIND) - - ! Initial displacement to skip at beginning of file - disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1) - - call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_DATA%view(i), & - 'native', mpi_info_int, ierr) - call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & - MPI_DOUBLE_PRECISION, status, ierr) - - end do - end if - - call s_mpi_barrier() - - call MPI_FILE_CLOSE(ifile, ierr) - - if (ib) then - - write (file_loc, '(A)') 'ib.dat' - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) - inquire (FILE=trim(file_loc), EXIST=file_exist) - - if (file_exist) then - - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - - disp = 0 - - call MPI_FILE_SET_VIEW(ifile, disp, MPI_INTEGER, MPI_IO_IB_DATA%view, & - 'native', mpi_info_int, ierr) - call MPI_FILE_READ(ifile, MPI_IO_IB_DATA%var%sf, data_size, & - MPI_INTEGER, status, ierr) - - else - call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') - end if - - end if - - else - call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') - end if - - end if - - if (ib) then - - do j = 1, num_ibs - if (patch_ib(j)%c > 0) then - - print *, "HERE Np", Np - - allocate (airfoil_grid_u(1:Np)) - allocate (airfoil_grid_l(1:Np)) - - write (file_loc, '(A)') 'airfoil_l.dat' - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) - inquire (FILE=trim(file_loc), EXIST=file_exist) - if (file_exist) then - - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - - ! Initial displacement to skip at beginning of file - disp = 0 - - call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_airfoil_IB_DATA%view(1), & - 'native', mpi_info_int, ierr) - call MPI_FILE_READ(ifile, MPI_IO_airfoil_IB_DATA%var(1:Np), 3*Np, & - MPI_DOUBLE_PRECISION, status, ierr) - - end if - - write (file_loc, '(A)') 'airfoil_u.dat' - file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) - inquire (FILE=trim(file_loc), EXIST=file_exist) - if (file_exist) then - - call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) - - ! Initial displacement to skip at beginning of file - disp = 0 - - call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_airfoil_IB_DATA%view(2), & - 'native', mpi_info_int, ierr) - call MPI_FILE_READ(ifile, MPI_IO_airfoil_IB_DATA%var(Np + 1:2*Np), 3*Np, & - MPI_DOUBLE_PRECISION, status, ierr) - end if - - do i = 1, Np - airfoil_grid_l(i)%x = MPI_IO_airfoil_IB_DATA%var(i)%x - airfoil_grid_l(i)%y = MPI_IO_airfoil_IB_DATA%var(i)%y - end do - - do i = 1, Np - airfoil_grid_u(i)%x = MPI_IO_airfoil_IB_DATA%var(Np + i)%x - airfoil_grid_u(i)%y = MPI_IO_airfoil_IB_DATA%var(Np + i)%y - end do - - end if - end do - end if - - deallocate (x_cb_glb, y_cb_glb, z_cb_glb) - -#endif - - end subroutine s_read_parallel_data_files - - !> The purpose of this subroutine is to open a new or pre- - !! existing run-time information file and append to it the - !! basic header information relevant to current simulation. - !! In general, this requires generating a table header for - !! those stability criteria which will be written at every - !! time-step. - subroutine s_open_run_time_information_file - - character(LEN=name_len) :: file_name = 'run_time.inf' !< - !! Name of the run-time information file - - character(LEN=path_len + name_len) :: file_path !< - !! Relative path to a file in the case directory - - character(LEN=8) :: file_date !< - !! Creation date of the run-time information file - - ! Opening the run-time information file - file_path = trim(case_dir)//'/'//trim(file_name) - - open (3, FILE=trim(file_path), & - FORM='formatted', & - STATUS='replace') - - write (3, '(A)') 'Description: Stability information at '// & - 'each time-step of the simulation. This' - write (3, '(13X,A)') 'data is composed of the inviscid '// & - 'Courant–Friedrichs–Lewy (ICFL)' - write (3, '(13X,A)') 'number, the viscous CFL (VCFL) number, '// & - 'the capillary CFL (CCFL)' - write (3, '(13X,A)') 'number and the cell Reynolds (Rc) '// & - 'number. Please note that only' - write (3, '(13X,A)') 'those stability conditions pertinent '// & - 'to the physics included in' - write (3, '(13X,A)') 'the current computation are displayed.' - - call date_and_time(DATE=file_date) - - write (3, '(A)') 'Date: '//file_date(5:6)//'/'// & - file_date(7:8)//'/'// & - file_date(3:4) - - write (3, '(A)') ''; write (3, '(A)') '' - - ! Generating table header for the stability criteria to be outputted - if (cfl_dt) then - if (any(Re_size > 0)) then - write (1, '(A)') '==== Time-steps ====== dt ===== Time ======= ICFL '// & - 'Max ==== VCFL Max ====== Rc Min =======' - else - write (1, '(A)') '=========== Time-steps ============== dt ===== Time '// & - '============== ICFL Max =============' - end if - else - if (any(Re_size > 0)) then - write (1, '(A)') '==== Time-steps ====== Time ======= ICFL '// & - 'Max ==== VCFL Max ====== Rc Min =======' - else - write (1, '(A)') '=========== Time-steps ============== Time '// & - '============== ICFL Max =============' - end if - end if - - end subroutine s_open_run_time_information_file - - !> The goal of the procedure is to output to the run-time - !! information file the stability criteria extrema in the - !! entire computational domain and at the given time-step. - !! Moreover, the subroutine is also in charge of tracking - !! these stability criteria extrema over all time-steps. - !! @param q_prim_vf Cell-average primitive variables - !! @param t_step Current time step - subroutine s_write_run_time_information(q_prim_vf, t_step) - - type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf - integer, intent(IN) :: t_step - - real(kind(0d0)), dimension(num_fluids) :: alpha_rho !< Cell-avg. partial density - real(kind(0d0)) :: rho !< Cell-avg. density - real(kind(0d0)), dimension(num_dims) :: vel !< Cell-avg. velocity - real(kind(0d0)) :: vel_sum !< Cell-avg. velocity sum - real(kind(0d0)) :: pres !< Cell-avg. pressure - real(kind(0d0)), dimension(num_fluids) :: alpha !< Cell-avg. volume fraction - real(kind(0d0)) :: gamma !< Cell-avg. sp. heat ratio - real(kind(0d0)) :: pi_inf !< Cell-avg. liquid stiffness function - real(kind(0d0)) :: qv !< Cell-avg. fluid reference energy - real(kind(0d0)) :: c !< Cell-avg. sound speed - real(kind(0d0)) :: E !< Cell-avg. energy - real(kind(0d0)) :: H !< Cell-avg. enthalpy - real(kind(0d0)), dimension(2) :: Re !< Cell-avg. Reynolds numbers - - ! ICFL, VCFL, CCFL and Rc stability criteria extrema for the current - ! time-step and located on both the local (loc) and the global (glb) - ! computational domains - - real(kind(0d0)) :: blkmod1, blkmod2 !< - !! Fluid bulk modulus for Woods mixture sound speed - - integer :: i, j, k, l, q !< Generic loop iterators - - integer :: Nfq - real(kind(0d0)) :: fltr_dtheta !< - !! Modified dtheta accounting for Fourier filtering in azimuthal direction. - - ! Computing Stability Criteria at Current Time-step ================ - !$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho, vel, alpha, Re, fltr_dtheta, Nfq) - do l = 0, p - do k = 0, n - do j = 0, m - - call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, j, k, l) - - ! Compute mixture sound speed - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, c) - - if (any(Re_size > 0)) then - call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf, vcfl_sf, Rc_sf) - else - call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf) - end if - - end do - end do - end do - ! END: Computing Stability Criteria at Current Time-step =========== - - ! Determining local stability criteria extrema at current time-step - -#ifdef CRAY_ACC_WAR - !$acc update host(icfl_sf) - - if (any(Re_size > 0)) then - !$acc update host(vcfl_sf, Rc_sf) - end if - - icfl_max_loc = maxval(icfl_sf) - - if (any(Re_size > 0)) then - vcfl_max_loc = maxval(vcfl_sf) - Rc_min_loc = minval(Rc_sf) - end if -#else - !$acc kernels - icfl_max_loc = maxval(icfl_sf) - !$acc end kernels - - if (any(Re_size > 0)) then - !$acc kernels - vcfl_max_loc = maxval(vcfl_sf) - Rc_min_loc = minval(Rc_sf) - !$acc end kernels - end if -#endif - - ! Determining global stability criteria extrema at current time-step - if (num_procs > 1) then - call s_mpi_reduce_stability_criteria_extrema(icfl_max_loc, & - vcfl_max_loc, & - ccfl_max_loc, & - Rc_min_loc, & - icfl_max_glb, & - vcfl_max_glb, & - ccfl_max_glb, & - Rc_min_glb) - else - icfl_max_glb = icfl_max_loc - if (any(Re_size > 0)) vcfl_max_glb = vcfl_max_loc - if (any(Re_size > 0)) Rc_min_glb = Rc_min_loc - end if - - ! Determining the stability criteria extrema over all the time-steps - if (icfl_max_glb > icfl_max) icfl_max = icfl_max_glb - - if (any(Re_size > 0)) then - if (vcfl_max_glb > vcfl_max) vcfl_max = vcfl_max_glb - if (Rc_min_glb < Rc_min) Rc_min = Rc_min_glb - end if - - ! Outputting global stability criteria extrema at current time-step - if (proc_rank == 0) then - if (any(Re_size > 0)) then - write (1, '(6X,I8,F10.6,6X,6X,F10.6,6X,F9.6,6X,F9.6,6X,F10.6)') & - t_step, dt, t_step*dt, icfl_max_glb, & - vcfl_max_glb, & - Rc_min_glb - else - write (1, '(13X,I8,14X,F10.6,14X,F10.6,13X,F9.6)') & - t_step, dt, t_step*dt, icfl_max_glb - end if - - if (icfl_max_glb /= icfl_max_glb) then - call s_mpi_abort('ICFL is NaN. Exiting ...') - elseif (icfl_max_glb > 1d0) then - print *, 'icfl', icfl_max_glb - call s_mpi_abort('ICFL is greater than 1.0. Exiting ...') - end if - - do i = chemxb, chemxe - !@:ASSERT(all(q_prim_vf(i)%sf(:,:,:) >= -1d0), "bad conc") - !@:ASSERT(all(q_prim_vf(i)%sf(:,:,:) <= 2d0), "bad conc") - end do - - if (vcfl_max_glb /= vcfl_max_glb) then - call s_mpi_abort('VCFL is NaN. Exiting ...') - elseif (vcfl_max_glb > 1d0) then - print *, 'vcfl', vcfl_max_glb - call s_mpi_abort('VCFL is greater than 1.0. Exiting ...') - end if - end if - - call s_mpi_barrier() - - end subroutine s_write_run_time_information - - !> The goal of this subroutine is to write to the run-time - !! information file basic footer information applicable to - !! the current computation and to close the file when done. - !! The footer contains the stability criteria extrema over - !! all of the time-steps and the simulation run-time. - subroutine s_close_run_time_information_file - - real(kind(0d0)) :: run_time !< Run-time of the simulation - ! Writing the footer of and closing the run-time information file - write (3, '(A)') '----------------------------------------'// & - '----------------------------------------' - write (3, '(A)') '' - - write (3, '(A,F9.6)') 'ICFL Max: ', icfl_max - if (any(Re_size > 0)) write (3, '(A,F9.6)') 'VCFL Max: ', vcfl_max - if (any(Re_size > 0)) write (3, '(A,F10.6)') 'Rc Min: ', Rc_min - - call cpu_time(run_time) - - write (3, '(A)') '' - write (3, '(A,I0,A)') 'Run-time: ', int(anint(run_time)), 's' - write (3, '(A)') '========================================'// & - '========================================' - close (3) - - end subroutine s_close_run_time_information_file - - subroutine s_finalize_sim_helpers_module() - - ! Deallocating the ICFL, VCFL, CCFL, and Rc stability criteria - @:DEALLOCATE_GLOBAL(icfl_sf) - if (any(Re_size > 0)) then - @:DEALLOCATE_GLOBAL(vcfl_sf, Rc_sf) - end if - - s_read_data_files => null() - s_write_data_files => null() - - end subroutine s_finalize_sim_helpers_module - -end module m_sim_helpers diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 05d9d2af0..004900eba 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -78,14 +78,15 @@ module m_start_up use m_surface_tension use m_body_forces - - use m_sim_helpers ! ========================================================================== implicit none private; public :: s_read_input_file, & s_check_input_file, & + s_read_data_files, & + s_read_serial_data_files, & + s_read_parallel_data_files, & s_populate_grid_variables_buffers, & s_initialize_internal_energy_equations, & s_initialize_modules, s_initialize_gpu_vars, & @@ -93,6 +94,25 @@ module m_start_up s_perform_time_step, s_save_data, & s_save_performance_metrics + abstract interface ! =================================================== + + !! @param q_cons_vf Conservative variables + subroutine s_read_abstract_data_files(q_cons_vf) + + import :: scalar_field, sys_size, pres_field + + type(scalar_field), & + dimension(sys_size), & + intent(inout) :: q_cons_vf + + end subroutine s_read_abstract_data_files + + end interface ! ======================================================== + + type(scalar_field), allocatable, dimension(:) :: grad_x_vf, grad_y_vf, grad_z_vf, norm_vf + + procedure(s_read_abstract_data_files), pointer :: s_read_data_files => null() + contains !> The purpose of this procedure is to first verify that an @@ -212,6 +232,596 @@ contains end subroutine s_check_input_file + !! initial condition and grid data files. The cell-average + !! conservative variables constitute the former, while the + !! cell-boundary locations in x-, y- and z-directions make + !! up the latter. This procedure also calculates the cell- + !! width distributions from the cell-boundary locations. + !! @param q_cons_vf Cell-averaged conservative variables + subroutine s_read_serial_data_files(q_cons_vf) + + type(scalar_field), dimension(sys_size), intent(INOUT) :: q_cons_vf + + character(LEN=path_len + 2*name_len) :: t_step_dir !< + !! Relative path to the starting time-step directory + + character(LEN=path_len + 3*name_len) :: file_path !< + !! Relative path to the grid and conservative variables data files + + logical :: file_exist !< + ! Logical used to check the existence of the data files + + integer :: i, r !< Generic loop iterator + + ! Confirming that the directory from which the initial condition and + ! the grid data files are to be read in exists and exiting otherwise + if (cfl_dt) then + write (t_step_dir, '(A,I0,A,I0)') & + trim(case_dir)//'/p_all/p', proc_rank, '/', n_start + else + write (t_step_dir, '(A,I0,A,I0)') & + trim(case_dir)//'/p_all/p', proc_rank, '/', t_step_start + end if + + file_path = trim(t_step_dir)//'/.' + call my_inquire(file_path, file_exist) + + if (file_exist .neqv. .true.) then + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + + ! Cell-boundary Locations in x-direction =========================== + file_path = trim(t_step_dir)//'/x_cb.dat' + + inquire (FILE=trim(file_path), EXIST=file_exist) + + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) x_cb(-1:m); close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + + dx(0:m) = x_cb(0:m) - x_cb(-1:m - 1) + x_cc(0:m) = x_cb(-1:m - 1) + dx(0:m)/2d0 + + if (ib) then + do i = 1, num_ibs + if (patch_ib(i)%c > 0) then + Np = int((patch_ib(i)%p*patch_ib(i)%c/dx(0))*20) + int(((patch_ib(i)%c - patch_ib(i)%p*patch_ib(i)%c)/dx(0))*20) + 1 + end if + end do + end if + ! ================================================================== + + ! Cell-boundary Locations in y-direction =========================== + if (n > 0) then + + file_path = trim(t_step_dir)//'/y_cb.dat' + + inquire (FILE=trim(file_path), EXIST=file_exist) + + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) y_cb(-1:n); close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + + dy(0:n) = y_cb(0:n) - y_cb(-1:n - 1) + y_cc(0:n) = y_cb(-1:n - 1) + dy(0:n)/2d0 + + end if + ! ================================================================== + + ! Cell-boundary Locations in z-direction =========================== + if (p > 0) then + + file_path = trim(t_step_dir)//'/z_cb.dat' + + inquire (FILE=trim(file_path), EXIST=file_exist) + + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) z_cb(-1:p); close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + + dz(0:p) = z_cb(0:p) - z_cb(-1:p - 1) + z_cc(0:p) = z_cb(-1:p - 1) + dz(0:p)/2d0 + + end if + ! ================================================================== + + do i = 1, sys_size + write (file_path, '(A,I0,A)') & + trim(t_step_dir)//'/q_cons_vf', i, '.dat' + inquire (FILE=trim(file_path), EXIST=file_exist) + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) q_cons_vf(i)%sf(0:m, 0:n, 0:p); close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + end do + + if ((bubbles .eqv. .true.) .or. (hypoelasticity .eqv. .true.)) then + ! Read pb and mv for non-polytropic qbmm + if (qbmm .and. .not. polytropic) then + do i = 1, nb + do r = 1, nnode + write (file_path, '(A,I0,A)') & + trim(t_step_dir)//'/pb', sys_size + (i - 1)*nnode + r, '.dat' + inquire (FILE=trim(file_path), EXIST=file_exist) + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) pb_ts(1)%sf(0:m, 0:n, 0:p, r, i); close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + end do + end do + do i = 1, nb + do r = 1, nnode + write (file_path, '(A,I0,A)') & + trim(t_step_dir)//'/mv', sys_size + (i - 1)*nnode + r, '.dat' + inquire (FILE=trim(file_path), EXIST=file_exist) + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) mv_ts(1)%sf(0:m, 0:n, 0:p, r, i); close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + end do + end do + end if + end if + ! ================================================================== + + ! Read IBM Data ==================================================== + + if (ib) then + do i = 1, num_ibs + write (file_path, '(A,I0,A)') & + trim(t_step_dir)//'/ib.dat' + inquire (FILE=trim(file_path), EXIST=file_exist) + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) ib_markers%sf(0:m, 0:n, 0:p); close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + + if (patch_ib(i)%c > 0) then + + print *, "HERE Np", Np + allocate (airfoil_grid_u(1:Np)) + allocate (airfoil_grid_l(1:Np)) + + write (file_path, '(A)') & + trim(t_step_dir)//'/airfoil_u.dat' + inquire (FILE=trim(file_path), EXIST=file_exist) + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) airfoil_grid_u; close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + + write (file_path, '(A)') & + trim(t_step_dir)//'/airfoil_l.dat' + inquire (FILE=trim(file_path), EXIST=file_exist) + if (file_exist) then + open (2, FILE=trim(file_path), & + FORM='unformatted', & + ACTION='read', & + STATUS='old') + read (2) airfoil_grid_l; close (2) + else + call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') + end if + end if + end do + + end if + + end subroutine s_read_serial_data_files + + !! @param q_cons_vf Conservative variables + subroutine s_read_parallel_data_files(q_cons_vf) + + type(scalar_field), & + dimension(sys_size), & + intent(INOUT) :: q_cons_vf + +#ifdef MFC_MPI + + real(kind(0d0)), allocatable, dimension(:) :: x_cb_glb, y_cb_glb, z_cb_glb + + 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) :: 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_start_string + + integer :: i, j + + allocate (x_cb_glb(-1:m_glb)) + allocate (y_cb_glb(-1:n_glb)) + allocate (z_cb_glb(-1:p_glb)) + + ! Read in cell boundary locations in x-direction + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//'x_cb.dat' + inquire (FILE=trim(file_loc), EXIST=file_exist) + + if (file_exist) then + data_size = m_glb + 2 + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + call MPI_FILE_READ(ifile, x_cb_glb, data_size, MPI_DOUBLE_PRECISION, status, ierr) + call MPI_FILE_CLOSE(ifile, ierr) + else + call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') + end if + + ! Assigning local cell boundary locations + x_cb(-1:m) = x_cb_glb((start_idx(1) - 1):(start_idx(1) + m)) + ! Computing the cell width distribution + dx(0:m) = x_cb(0:m) - x_cb(-1:m - 1) + ! Computing the cell center locations + x_cc(0:m) = x_cb(-1:m - 1) + dx(0:m)/2d0 + + if (ib) then + do i = 1, num_ibs + if (patch_ib(i)%c > 0) then + Np = int((patch_ib(i)%p*patch_ib(i)%c/dx(0))*20) + int(((patch_ib(i)%c - patch_ib(i)%p*patch_ib(i)%c)/dx(0))*20) + 1 + allocate (MPI_IO_airfoil_IB_DATA%var(1:2*Np)) + print *, "HERE Np", Np + end if + end do + end if + + if (n > 0) then + ! Read in cell boundary locations in y-direction + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//'y_cb.dat' + inquire (FILE=trim(file_loc), EXIST=file_exist) + + if (file_exist) then + data_size = n_glb + 2 + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + call MPI_FILE_READ(ifile, y_cb_glb, data_size, MPI_DOUBLE_PRECISION, status, ierr) + call MPI_FILE_CLOSE(ifile, ierr) + else + call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') + end if + + ! Assigning local cell boundary locations + y_cb(-1:n) = y_cb_glb((start_idx(2) - 1):(start_idx(2) + n)) + ! Computing the cell width distribution + dy(0:n) = y_cb(0:n) - y_cb(-1:n - 1) + ! Computing the cell center locations + y_cc(0:n) = y_cb(-1:n - 1) + dy(0:n)/2d0 + + if (p > 0) then + ! Read in cell boundary locations in z-direction + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//'z_cb.dat' + inquire (FILE=trim(file_loc), EXIST=file_exist) + + if (file_exist) then + data_size = p_glb + 2 + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + call MPI_FILE_READ(ifile, z_cb_glb, data_size, MPI_DOUBLE_PRECISION, status, ierr) + call MPI_FILE_CLOSE(ifile, ierr) + else + call s_mpi_abort('File '//trim(file_loc)//'is missing. Exiting...') + end if + + ! Assigning local cell boundary locations + z_cb(-1:p) = z_cb_glb((start_idx(3) - 1):(start_idx(3) + p)) + ! Computing the cell width distribution + dz(0:p) = z_cb(0:p) - z_cb(-1:p - 1) + ! Computing the cell center locations + z_cc(0:p) = z_cb(-1:p - 1) + dz(0:p)/2d0 + + end if + end if + + if (file_per_process) then + if (cfl_dt) then + call s_int_to_str(n_start, t_step_start_string) + write (file_loc, '(I0,A1,I7.7,A)') n_start, '_', proc_rank, '.dat' + else + call s_int_to_str(t_step_start, t_step_start_string) + write (file_loc, '(I0,A1,I7.7,A)') t_step_start, '_', proc_rank, '.dat' + end if + + file_loc = trim(case_dir)//'/restart_data/lustre_'//trim(t_step_start_string)//trim(mpiiofs)//trim(file_loc) + inquire (FILE=trim(file_loc), EXIST=file_exist) + + if (file_exist) then + call MPI_FILE_OPEN(MPI_COMM_SELF, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + + ! Initialize MPI data I/O + + if (ib) then + call s_initialize_mpi_data(q_cons_vf, ib_markers) + else + call s_initialize_mpi_data(q_cons_vf) + end if + + ! Size of local arrays + data_size = (m + 1)*(n + 1)*(p + 1) + + ! Resize some integers so MPI can read even the biggest file + m_MOK = int(m_glb + 1, MPI_OFFSET_KIND) + n_MOK = int(n_glb + 1, MPI_OFFSET_KIND) + p_MOK = int(p_glb + 1, MPI_OFFSET_KIND) + WP_MOK = int(8d0, MPI_OFFSET_KIND) + MOK = int(1d0, MPI_OFFSET_KIND) + str_MOK = int(name_len, MPI_OFFSET_KIND) + NVARS_MOK = int(sys_size, MPI_OFFSET_KIND) + + ! Read the data for each variable + if (bubbles .or. hypoelasticity) then + + do i = 1, sys_size!adv_idx%end + var_MOK = int(i, MPI_OFFSET_KIND) + + call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & + MPI_DOUBLE_PRECISION, status, ierr) + end do + !Read pb and mv for non-polytropic qbmm + if (qbmm .and. .not. polytropic) then + do i = sys_size + 1, sys_size + 2*nb*nnode + var_MOK = int(i, MPI_OFFSET_KIND) + + call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & + MPI_DOUBLE_PRECISION, status, ierr) + end do + end if + else + do i = 1, adv_idx%end + var_MOK = int(i, MPI_OFFSET_KIND) + + call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & + MPI_DOUBLE_PRECISION, status, ierr) + end do + end if + + call s_mpi_barrier() + + call MPI_FILE_CLOSE(ifile, ierr) + + if (ib) then + + write (file_loc, '(A)') 'ib.dat' + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) + inquire (FILE=trim(file_loc), EXIST=file_exist) + + if (file_exist) then + + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + + disp = 0 + + call MPI_FILE_SET_VIEW(ifile, disp, MPI_INTEGER, MPI_IO_IB_DATA%view, & + 'native', mpi_info_int, ierr) + call MPI_FILE_READ(ifile, MPI_IO_IB_DATA%var%sf, data_size, & + MPI_INTEGER, status, ierr) + + else + call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') + end if + + end if + + else + call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') + end if + else + ! Open the file to read conservative variables + if (cfl_dt) then + write (file_loc, '(I0,A)') n_start, '.dat' + else + write (file_loc, '(I0,A)') t_step_start, '.dat' + end if + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) + inquire (FILE=trim(file_loc), EXIST=file_exist) + + if (file_exist) then + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + + ! Initialize MPI data I/O + + if (ib) then + call s_initialize_mpi_data(q_cons_vf, ib_markers) + else + + call s_initialize_mpi_data(q_cons_vf) + + end if + + + ! Size of local arrays + data_size = (m + 1)*(n + 1)*(p + 1) + + ! Resize some integers so MPI can read even the biggest file + m_MOK = int(m_glb + 1, MPI_OFFSET_KIND) + n_MOK = int(n_glb + 1, MPI_OFFSET_KIND) + p_MOK = int(p_glb + 1, MPI_OFFSET_KIND) + WP_MOK = int(8d0, MPI_OFFSET_KIND) + MOK = int(1d0, MPI_OFFSET_KIND) + str_MOK = int(name_len, MPI_OFFSET_KIND) + NVARS_MOK = int(sys_size, MPI_OFFSET_KIND) + + ! Read the data for each variable + if (bubbles .or. hypoelasticity) then + + do i = 1, sys_size!adv_idx%end + var_MOK = int(i, MPI_OFFSET_KIND) + ! Initial displacement to skip at beginning of file + disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1) + + call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_DATA%view(i), & + 'native', mpi_info_int, ierr) + call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & + MPI_DOUBLE_PRECISION, status, ierr) + end do + !Read pb and mv for non-polytropic qbmm + if (qbmm .and. .not. polytropic) then + do i = sys_size + 1, sys_size + 2*nb*nnode + var_MOK = int(i, MPI_OFFSET_KIND) + ! Initial displacement to skip at beginning of file + disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1) + + call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_DATA%view(i), & + 'native', mpi_info_int, ierr) + call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & + MPI_DOUBLE_PRECISION, status, ierr) + end do + end if + else + do i = 1, sys_size + var_MOK = int(i, MPI_OFFSET_KIND) + + ! Initial displacement to skip at beginning of file + disp = m_MOK*max(MOK, n_MOK)*max(MOK, p_MOK)*WP_MOK*(var_MOK - 1) + + call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_DATA%view(i), & + 'native', mpi_info_int, ierr) + call MPI_FILE_READ(ifile, MPI_IO_DATA%var(i)%sf, data_size, & + MPI_DOUBLE_PRECISION, status, ierr) + + end do + end if + + call s_mpi_barrier() + + call MPI_FILE_CLOSE(ifile, ierr) + + if (ib) then + + write (file_loc, '(A)') 'ib.dat' + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) + inquire (FILE=trim(file_loc), EXIST=file_exist) + + if (file_exist) then + + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + + disp = 0 + + call MPI_FILE_SET_VIEW(ifile, disp, MPI_INTEGER, MPI_IO_IB_DATA%view, & + 'native', mpi_info_int, ierr) + call MPI_FILE_READ(ifile, MPI_IO_IB_DATA%var%sf, data_size, & + MPI_INTEGER, status, ierr) + + else + call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') + end if + + end if + + else + call s_mpi_abort('File '//trim(file_loc)//' is missing. Exiting...') + end if + + end if + + if (ib) then + + do j = 1, num_ibs + if (patch_ib(j)%c > 0) then + + print *, "HERE Np", Np + + allocate (airfoil_grid_u(1:Np)) + allocate (airfoil_grid_l(1:Np)) + + write (file_loc, '(A)') 'airfoil_l.dat' + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) + inquire (FILE=trim(file_loc), EXIST=file_exist) + if (file_exist) then + + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + + ! Initial displacement to skip at beginning of file + disp = 0 + + call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_airfoil_IB_DATA%view(1), & + 'native', mpi_info_int, ierr) + call MPI_FILE_READ(ifile, MPI_IO_airfoil_IB_DATA%var(1:Np), 3*Np, & + MPI_DOUBLE_PRECISION, status, ierr) + + end if + + write (file_loc, '(A)') 'airfoil_u.dat' + file_loc = trim(case_dir)//'/restart_data'//trim(mpiiofs)//trim(file_loc) + inquire (FILE=trim(file_loc), EXIST=file_exist) + if (file_exist) then + + call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, MPI_MODE_RDONLY, mpi_info_int, ifile, ierr) + + ! Initial displacement to skip at beginning of file + disp = 0 + + call MPI_FILE_SET_VIEW(ifile, disp, MPI_DOUBLE_PRECISION, MPI_IO_airfoil_IB_DATA%view(2), & + 'native', mpi_info_int, ierr) + call MPI_FILE_READ(ifile, MPI_IO_airfoil_IB_DATA%var(Np + 1:2*Np), 3*Np, & + MPI_DOUBLE_PRECISION, status, ierr) + end if + + do i = 1, Np + airfoil_grid_l(i)%x = MPI_IO_airfoil_IB_DATA%var(i)%x + airfoil_grid_l(i)%y = MPI_IO_airfoil_IB_DATA%var(i)%y + end do + + do i = 1, Np + airfoil_grid_u(i)%x = MPI_IO_airfoil_IB_DATA%var(Np + i)%x + airfoil_grid_u(i)%y = MPI_IO_airfoil_IB_DATA%var(Np + i)%y + end do + + end if + end do + end if + + deallocate (x_cb_glb, y_cb_glb, z_cb_glb) + +#endif + + end subroutine s_read_parallel_data_files + !> The purpose of this subroutine is to populate the buffers !! of the grid variables, which are constituted of the cell- !! boundary locations and cell-width distributions, based on @@ -696,8 +1306,6 @@ contains pref = 1d0 end if - call s_initialize_sim_helpers_module() - #if defined(MFC_OpenACC) && defined(MFC_MEMORY_DUMP) call acc_present_dump() #endif @@ -743,8 +1351,17 @@ contains call acc_present_dump() #endif + ! Associate pointers for serial or parallel I/O + if (parallel_io .neqv. .true.) then + s_read_data_files => s_read_serial_data_files + s_write_data_files => s_write_serial_data_files + else + s_read_data_files => s_read_parallel_data_files + s_write_data_files => s_write_parallel_data_files + end if + ! Reading in the user provided initial condition and grid data - call s_read_data_files(q_cons_ts(1)%vf, ib_markers) + call s_read_data_files(q_cons_ts(1)%vf) if (model_eqns == 3) call s_initialize_internal_energy_equations(q_cons_ts(1)%vf) if (ib) call s_ibm_setup() @@ -880,7 +1497,6 @@ contains s_read_data_files => null() s_write_data_files => null() - call s_finalize_sim_helpers_module() call s_finalize_time_steppers_module() call s_finalize_derived_variables_module() call s_finalize_data_output_module() diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index f7a90e6da..e7327190e 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -31,6 +31,8 @@ module m_time_steppers use m_helper + use m_sim_helpers + use m_fftw use m_nvtx @@ -1028,6 +1030,95 @@ contains end subroutine s_apply_bodyforces + subroutine s_comprehensive_debug(q_cons_vf, q_prim_vf, t_step, stage) + + type(scalar_field), dimension(sys_size) :: q_cons_vf, q_prim_vf + integer, intent(in) :: t_step, stage + integer :: j, k, l, i + integer errors + logical :: exists + + character(LEN=name_len) :: file_name = 'comp_debug.txt' + character(LEN=path_len + name_len) :: file_path + character(100) :: str_format + + ! Opening the run-time information file + file_path = trim(case_dir)//'/'//trim(file_name) + + str_format = "(I9, A, I3, A, I4, I4, I4, A, I2, A, I5, A, I5, I5, I5)" + + open (12, FILE=trim(file_path), & + STATUS='replace') + + errors = 0 + + ! Check all variables for NaNs + do i = 1, sys_size + !$acc update host(q_cons_vf(i)%sf) + do l = 0, p + do k = 0, n + do j = 0, m + if (ieee_is_nan(q_cons_vf(i)%sf(j, k, l))) then + write(12, str_format) t_step, " NaN(s) in conservative variables after RK stage ", & + stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & + " (m, n, p)", m, n, p + errors = errors + 1 + end if + end do + end do + end do + end do + + ! Check for invalid volume fractions + do i = advxb, advxe + do l = 0, p + do k = 0, n + do j = 0, m + if (q_cons_vf(i)%sf(j, k, l) < 0d0) then + write(12, str_format) t_step, " Volume fraction < 0 after RK stage ", & + stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & + " (m, n, p)", m, n, p + errors = errors + 1 + elseif (q_cons_vf(i)%sf(j, k, l) > 1d0 + verysmall) then + write(12, str_format) t_step, " Volume fraction > 1 after RK stage ", & + stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & + " (m, n, p)", m, n, p + errors = errors + 1 + end if + end do + end do + end do + end do + + ! Check for invalid densities + do i = contxb, contxe + do l = 0, p + do k = 0, n + do j = 0, m + if (q_cons_vf(advxb + i -1)%sf(j, k, l) < 0d0 .and. q_cons_vf(i)%sf(j, k, l) < 0d0 .or. & + q_cons_vf(advxb + i -1)%sf(j, k, l) > 0d0 .and. q_cons_Vf(i)%sf(j, k, l) < 0d0) then + print*, q_cons_vf(advxb + i - 1)%sf(j, k, l), q_cons_vf(i)%sf(j, k, l) + write(12, str_format) t_step, " Density is negative after RK stage ", & + stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & + " (m, n, p)", m, n, p + errors = errors + 1 + end if + end do + end do + end do + end do + + if (errors /= 0) then + close(12) + call s_write_data_files(q_cons_vf, q_prim_vf, t_step) + call s_mpi_abort("Errors found in conservative variables") + endif + + write(12, "(I3)") -1 + close(12) + + end subroutine s_comprehensive_debug + !> This subroutine saves the temporary q_prim_vf vector !! into the q_prim_ts vector that is then used in p_main !! @param t_step current time-step From dfe522cf7a9dfb8f64b74debb692c09a50a0241b Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Wed, 16 Oct 2024 14:20:41 -0400 Subject: [PATCH 07/10] update things --- src/simulation/m_sim_helpers.f90 | 83 +++++++++++++++++++++++++++++- src/simulation/m_time_steppers.fpp | 77 +++------------------------ 2 files changed, 88 insertions(+), 72 deletions(-) diff --git a/src/simulation/m_sim_helpers.f90 b/src/simulation/m_sim_helpers.f90 index be95197f6..a0e425d06 100644 --- a/src/simulation/m_sim_helpers.f90 +++ b/src/simulation/m_sim_helpers.f90 @@ -14,7 +14,8 @@ module m_sim_helpers private; public :: s_compute_enthalpy, & s_compute_stability_from_dt, & - s_compute_dt_from_cfl + s_compute_dt_from_cfl, & + s_check_cells contains @@ -247,4 +248,84 @@ subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l) end subroutine s_compute_dt_from_cfl + subroutine s_check_cells(q_cons_Vf, q_prim_Vf, t_step, stage, errors) + + type(scalar_field), dimension(sys_size) :: q_cons_vf, q_prim_vf + integer, intent(in) :: t_step, stage + integer :: j, k, l, i + integer errors + logical :: exists + + character(LEN=name_len) :: file_name = 'comp_debug.txt' + character(LEN=path_len + name_len) :: file_path + character(100) :: str_format + + ! Opening the run-time information file + file_path = trim(case_dir)//'/'//trim(file_name) + + str_format = "(I9, A, I3, A, I4, I4, I4, A, I2, A, I5, A, I5, I5, I5)" + + open (12, FILE=trim(file_path), & + STATUS='replace') + + errors = 0 + + ! Check all variables for NaNs + do i = 1, sys_size + !$acc update host(q_cons_vf(i)%sf) + do l = 0, p + do k = 0, n + do j = 0, m + if (ieee_is_nan(q_cons_vf(i)%sf(j, k, l))) then + write(12, str_format) t_step, " NaN(s) in conservative variables after RK stage ", & + stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & + " (m, n, p)", m, n, p + errors = errors + 1 + end if + end do + end do + end do + end do + + ! Check for invalid volume fractions + do i = advxb, advxe + do l = 0, p + do k = 0, n + do j = 0, m + if (q_cons_vf(i)%sf(j, k, l) < 0d0) then + write(12, str_format) t_step, " Volume fraction < 0 after RK stage ", & + stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & + " (m, n, p)", m, n, p + errors = errors + 1 + elseif (q_cons_vf(i)%sf(j, k, l) > 1d0 + verysmall) then + write(12, str_format) t_step, " Volume fraction > 1 after RK stage ", & + stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & + " (m, n, p)", m, n, p + errors = errors + 1 + end if + end do + end do + end do + end do + + ! Check for invalid densities + do i = contxb, contxe + do l = 0, p + do k = 0, n + do j = 0, m + if (q_cons_vf(advxb + i -1)%sf(j, k, l) < 0d0 .and. q_cons_vf(i)%sf(j, k, l) < 0d0 .or. & + q_cons_vf(advxb + i -1)%sf(j, k, l) > 0d0 .and. q_cons_Vf(i)%sf(j, k, l) < 0d0) then + print*, q_cons_vf(advxb + i - 1)%sf(j, k, l), q_cons_vf(i)%sf(j, k, l) + write(12, str_format) t_step, " Density is negative after RK stage ", & + stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & + " (m, n, p)", m, n, p + errors = errors + 1 + end if + end do + end do + end do + end do + + end subroutine + end module m_sim_helpers diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index e7327190e..c2a430bfd 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -866,6 +866,10 @@ contains call s_pressure_relaxation_procedure(q_cons_ts(1)%vf) end if + q_cons_ts(1)%vf(advxb)%sf(10, 10, 0) = -1 + q_cons_ts(1)%vf(contxe)%sf(10, 10, 0) = -1 + q_Cons_ts(1)%vf(advxe)%sf(5, 5, 0) = 1.1 + if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, q_prim_vf, t_step, 3) if (adv_n) call s_comp_alpha_from_n(q_cons_ts(1)%vf) @@ -1034,79 +1038,10 @@ contains type(scalar_field), dimension(sys_size) :: q_cons_vf, q_prim_vf integer, intent(in) :: t_step, stage - integer :: j, k, l, i - integer errors - logical :: exists - - character(LEN=name_len) :: file_name = 'comp_debug.txt' - character(LEN=path_len + name_len) :: file_path - character(100) :: str_format - - ! Opening the run-time information file - file_path = trim(case_dir)//'/'//trim(file_name) - - str_format = "(I9, A, I3, A, I4, I4, I4, A, I2, A, I5, A, I5, I5, I5)" - open (12, FILE=trim(file_path), & - STATUS='replace') + integer :: errors - errors = 0 - - ! Check all variables for NaNs - do i = 1, sys_size - !$acc update host(q_cons_vf(i)%sf) - do l = 0, p - do k = 0, n - do j = 0, m - if (ieee_is_nan(q_cons_vf(i)%sf(j, k, l))) then - write(12, str_format) t_step, " NaN(s) in conservative variables after RK stage ", & - stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & - " (m, n, p)", m, n, p - errors = errors + 1 - end if - end do - end do - end do - end do - - ! Check for invalid volume fractions - do i = advxb, advxe - do l = 0, p - do k = 0, n - do j = 0, m - if (q_cons_vf(i)%sf(j, k, l) < 0d0) then - write(12, str_format) t_step, " Volume fraction < 0 after RK stage ", & - stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & - " (m, n, p)", m, n, p - errors = errors + 1 - elseif (q_cons_vf(i)%sf(j, k, l) > 1d0 + verysmall) then - write(12, str_format) t_step, " Volume fraction > 1 after RK stage ", & - stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & - " (m, n, p)", m, n, p - errors = errors + 1 - end if - end do - end do - end do - end do - - ! Check for invalid densities - do i = contxb, contxe - do l = 0, p - do k = 0, n - do j = 0, m - if (q_cons_vf(advxb + i -1)%sf(j, k, l) < 0d0 .and. q_cons_vf(i)%sf(j, k, l) < 0d0 .or. & - q_cons_vf(advxb + i -1)%sf(j, k, l) > 0d0 .and. q_cons_Vf(i)%sf(j, k, l) < 0d0) then - print*, q_cons_vf(advxb + i - 1)%sf(j, k, l), q_cons_vf(i)%sf(j, k, l) - write(12, str_format) t_step, " Density is negative after RK stage ", & - stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & - " (m, n, p)", m, n, p - errors = errors + 1 - end if - end do - end do - end do - end do + call s_check_cells(q_cons_vf, q_prim_vf, t_step, stage, errors) if (errors /= 0) then close(12) From 817d338e6bdb00ea7462b707102c30716bd39f5d Mon Sep 17 00:00:00 2001 From: wilfonba Date: Wed, 16 Oct 2024 16:43:34 -0400 Subject: [PATCH 08/10] cleaning up --- examples/3D_performance_test/comp_debug.txt | 3 +++ src/post_process/m_start_up.f90 | 16 ++++++------ src/simulation/m_global_parameters.fpp | 8 +++--- src/simulation/m_sim_helpers.f90 | 20 +++++++-------- src/simulation/m_time_steppers.fpp | 28 +++++++++------------ 5 files changed, 37 insertions(+), 38 deletions(-) create mode 100644 examples/3D_performance_test/comp_debug.txt diff --git a/examples/3D_performance_test/comp_debug.txt b/examples/3D_performance_test/comp_debug.txt new file mode 100644 index 000000000..93cd8fd53 --- /dev/null +++ b/examples/3D_performance_test/comp_debug.txt @@ -0,0 +1,3 @@ + 0 Volume fraction < 0 after RK stage 2 at (j,k,l) 10 10 0 equation 7 proc 0 (m, n, p) 200 200 200 + 0 Volume fraction > 1 after RK stage 2 at (j,k,l) 5 5 0 equation 8 proc 0 (m, n, p) 200 200 200 + 0 Density is negative after RK stage 2 at (j,k,l) 10 10 0 equation 2 proc 0 (m, n, p) 200 200 200 diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90 index 3de39b661..9709e9387 100644 --- a/src/post_process/m_start_up.f90 +++ b/src/post_process/m_start_up.f90 @@ -198,15 +198,15 @@ subroutine s_perform_comprehensive_debug(varname, pres, c, H) ! Opening the run-time information file file_path = trim(case_dir)//'/'//trim(file_name) - inquire(file=file_path, exist=exists) + inquire (file=file_path, exist=exists) if (exists) then - open(12, file=file_path) - ! Read the file line by line - do - read(12, '(I9)', iostat=ios) t_fail - if (ios /= 0) exit ! Exit loop on error or end-of-file - end do - if (t_fail == -1) return + open (12, file=file_path) + ! Read the file line by line + do + read (12, '(I9)', iostat=ios) t_fail + if (ios /= 0) exit ! Exit loop on error or end-of-file + end do + if (t_fail == -1) return else return end if diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index dde62a08f..e4aa00037 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -1120,12 +1120,12 @@ contains @:ALLOCATE_GLOBAL(x_cc(-buff_size:m + buff_size)) @:ALLOCATE_GLOBAL(dx(-buff_size:m + buff_size)) - if (n == 0) return; + if (n == 0) return; @:ALLOCATE_GLOBAL(y_cb(-1 - buff_size:n + buff_size)) @:ALLOCATE_GLOBAL(y_cc(-buff_size:n + buff_size)) @:ALLOCATE_GLOBAL(dy(-buff_size:n + buff_size)) - if (p == 0) return; + if (p == 0) return; @:ALLOCATE_GLOBAL(z_cb(-1 - buff_size:p + buff_size)) @:ALLOCATE_GLOBAL(z_cc(-buff_size:p + buff_size)) @:ALLOCATE_GLOBAL(dz(-buff_size:p + buff_size)) @@ -1191,10 +1191,10 @@ contains ! Deallocating grid variables for the x-, y- and z-directions @:DEALLOCATE_GLOBAL(x_cb, x_cc, dx) - if (n == 0) return; + if (n == 0) return; @:DEALLOCATE_GLOBAL(y_cb, y_cc, dy) - if (p == 0) return; + if (p == 0) return; @:DEALLOCATE_GLOBAL(z_cb, z_cc, dz) end subroutine s_finalize_global_parameters_module diff --git a/src/simulation/m_sim_helpers.f90 b/src/simulation/m_sim_helpers.f90 index a0e425d06..cfbc19985 100644 --- a/src/simulation/m_sim_helpers.f90 +++ b/src/simulation/m_sim_helpers.f90 @@ -248,9 +248,9 @@ subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l) end subroutine s_compute_dt_from_cfl - subroutine s_check_cells(q_cons_Vf, q_prim_Vf, t_step, stage, errors) + subroutine s_check_cells(q_cons_Vf, t_step, stage, errors) - type(scalar_field), dimension(sys_size) :: q_cons_vf, q_prim_vf + type(scalar_field), dimension(sys_size) :: q_cons_vf integer, intent(in) :: t_step, stage integer :: j, k, l, i integer errors @@ -266,7 +266,7 @@ subroutine s_check_cells(q_cons_Vf, q_prim_Vf, t_step, stage, errors) str_format = "(I9, A, I3, A, I4, I4, I4, A, I2, A, I5, A, I5, I5, I5)" open (12, FILE=trim(file_path), & - STATUS='replace') + STATUS='replace') errors = 0 @@ -277,7 +277,7 @@ subroutine s_check_cells(q_cons_Vf, q_prim_Vf, t_step, stage, errors) do k = 0, n do j = 0, m if (ieee_is_nan(q_cons_vf(i)%sf(j, k, l))) then - write(12, str_format) t_step, " NaN(s) in conservative variables after RK stage ", & + write (12, str_format) t_step, " NaN(s) in conservative variables after RK stage ", & stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & " (m, n, p)", m, n, p errors = errors + 1 @@ -293,12 +293,12 @@ subroutine s_check_cells(q_cons_Vf, q_prim_Vf, t_step, stage, errors) do k = 0, n do j = 0, m if (q_cons_vf(i)%sf(j, k, l) < 0d0) then - write(12, str_format) t_step, " Volume fraction < 0 after RK stage ", & + write (12, str_format) t_step, " Volume fraction < 0 after RK stage ", & stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & " (m, n, p)", m, n, p errors = errors + 1 elseif (q_cons_vf(i)%sf(j, k, l) > 1d0 + verysmall) then - write(12, str_format) t_step, " Volume fraction > 1 after RK stage ", & + write (12, str_format) t_step, " Volume fraction > 1 after RK stage ", & stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & " (m, n, p)", m, n, p errors = errors + 1 @@ -313,10 +313,10 @@ subroutine s_check_cells(q_cons_Vf, q_prim_Vf, t_step, stage, errors) do l = 0, p do k = 0, n do j = 0, m - if (q_cons_vf(advxb + i -1)%sf(j, k, l) < 0d0 .and. q_cons_vf(i)%sf(j, k, l) < 0d0 .or. & - q_cons_vf(advxb + i -1)%sf(j, k, l) > 0d0 .and. q_cons_Vf(i)%sf(j, k, l) < 0d0) then - print*, q_cons_vf(advxb + i - 1)%sf(j, k, l), q_cons_vf(i)%sf(j, k, l) - write(12, str_format) t_step, " Density is negative after RK stage ", & + if (q_cons_vf(advxb + i - 1)%sf(j, k, l) < 0d0 .and. q_cons_vf(i)%sf(j, k, l) < 0d0 .or. & + q_cons_vf(advxb + i - 1)%sf(j, k, l) > 0d0 .and. q_cons_Vf(i)%sf(j, k, l) < 0d0) then + print *, q_cons_vf(advxb + i - 1)%sf(j, k, l), q_cons_vf(i)%sf(j, k, l) + write (12, str_format) t_step, " Density is negative after RK stage ", & stage, " at (j,k,l) ", j, k, l, " equation", i, " proc", proc_rank, & " (m, n, p)", m, n, p errors = errors + 1 diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index c2a430bfd..c853a1a9b 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -415,7 +415,7 @@ contains if (model_eqns == 3) call s_pressure_relaxation_procedure(q_cons_ts(1)%vf) - if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, q_prim_vf, t_step, 1) + if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, t_step, 1) if (adv_n) call s_comp_alpha_from_n(q_cons_ts(1)%vf) @@ -524,7 +524,7 @@ contains call s_pressure_relaxation_procedure(q_cons_ts(2)%vf) end if - if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, q_prim_vf, t_step, 1) + if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, t_step, 1) if (adv_n) call s_comp_alpha_from_n(q_cons_ts(2)%vf) @@ -601,7 +601,7 @@ contains call s_pressure_relaxation_procedure(q_cons_ts(1)%vf) end if - if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, q_prim_vf, t_step, 2) + if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, t_step, 2) if (adv_n) call s_comp_alpha_from_n(q_cons_ts(1)%vf) @@ -713,7 +713,7 @@ contains call s_pressure_relaxation_procedure(q_cons_ts(2)%vf) end if - if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, q_prim_vf, t_step, 1) + if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, t_step, 1) if (adv_n) call s_comp_alpha_from_n(q_cons_ts(2)%vf) @@ -790,7 +790,7 @@ contains call s_pressure_relaxation_procedure(q_cons_ts(2)%vf) end if - if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, q_prim_vf, t_step, 2) + if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, t_step, 2) if (adv_n) call s_comp_alpha_from_n(q_cons_ts(2)%vf) @@ -866,11 +866,7 @@ contains call s_pressure_relaxation_procedure(q_cons_ts(1)%vf) end if - q_cons_ts(1)%vf(advxb)%sf(10, 10, 0) = -1 - q_cons_ts(1)%vf(contxe)%sf(10, 10, 0) = -1 - q_Cons_ts(1)%vf(advxe)%sf(5, 5, 0) = 1.1 - - if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, q_prim_vf, t_step, 3) + if (comp_debug) call s_comprehensive_debug(q_cons_ts(1)%vf, t_step, 3) if (adv_n) call s_comp_alpha_from_n(q_cons_ts(1)%vf) @@ -1034,9 +1030,9 @@ contains end subroutine s_apply_bodyforces - subroutine s_comprehensive_debug(q_cons_vf, q_prim_vf, t_step, stage) + subroutine s_comprehensive_debug(q_cons_vf, t_step, stage) - type(scalar_field), dimension(sys_size) :: q_cons_vf, q_prim_vf + type(scalar_field), dimension(sys_size) :: q_cons_vf integer, intent(in) :: t_step, stage integer :: errors @@ -1044,13 +1040,13 @@ contains call s_check_cells(q_cons_vf, q_prim_vf, t_step, stage, errors) if (errors /= 0) then - close(12) + close (12) call s_write_data_files(q_cons_vf, q_prim_vf, t_step) call s_mpi_abort("Errors found in conservative variables") - endif + end if - write(12, "(I3)") -1 - close(12) + write (12, "(I3)") - 1 + close (12) end subroutine s_comprehensive_debug From f885194bd0ed6095cfd3e502d5392bf7dd8c013a Mon Sep 17 00:00:00 2001 From: wilfonba Date: Wed, 16 Oct 2024 17:26:46 -0400 Subject: [PATCH 09/10] cleaning up --- examples/3D_performance_test/comp_debug.txt | 4 +--- src/simulation/m_time_steppers.fpp | 2 +- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/examples/3D_performance_test/comp_debug.txt b/examples/3D_performance_test/comp_debug.txt index 93cd8fd53..af694311d 100644 --- a/examples/3D_performance_test/comp_debug.txt +++ b/examples/3D_performance_test/comp_debug.txt @@ -1,3 +1 @@ - 0 Volume fraction < 0 after RK stage 2 at (j,k,l) 10 10 0 equation 7 proc 0 (m, n, p) 200 200 200 - 0 Volume fraction > 1 after RK stage 2 at (j,k,l) 5 5 0 equation 8 proc 0 (m, n, p) 200 200 200 - 0 Density is negative after RK stage 2 at (j,k,l) 10 10 0 equation 2 proc 0 (m, n, p) 200 200 200 + -1 diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index c853a1a9b..4e9fabfc9 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -1037,7 +1037,7 @@ contains integer :: errors - call s_check_cells(q_cons_vf, q_prim_vf, t_step, stage, errors) + call s_check_cells(q_cons_vf, t_step, stage, errors) if (errors /= 0) then close (12) From 7c4842888b44de12c406d3e23563047317cafae3 Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Wed, 16 Oct 2024 17:54:27 -0400 Subject: [PATCH 10/10] fix cleanliness --- src/simulation/m_time_steppers.fpp | 13 +------------ 1 file changed, 1 insertion(+), 12 deletions(-) diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 4e9fabfc9..9ca06fb07 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -325,11 +325,6 @@ contains real(kind(0d0)), intent(inout) :: time_avg integer :: i, j, k, l, q!< Generic loop iterator - real(kind(0d0)) :: nR3bar - real(kind(0d0)) :: e_mix - - real(kind(0d0)) :: T - real(kind(0d0)), dimension(num_species) :: Ys ! Stage 1 of 1 ===================================================== @@ -442,7 +437,6 @@ contains integer :: i, j, k, l, q!< Generic loop iterator real(kind(0d0)) :: start, finish - real(kind(0d0)) :: nR3bar ! Stage 1 of 2 ===================================================== @@ -628,9 +622,7 @@ contains real(kind(0d0)), intent(INOUT) :: time_avg integer :: i, j, k, l, q !< Generic loop iterator - real(kind(0d0)) :: ts_error, denom, error_fraction, time_step_factor !< Generic loop iterator real(kind(0d0)) :: start, finish - real(kind(0d0)) :: nR3bar ! Stage 1 of 3 ===================================================== @@ -897,7 +889,6 @@ contains integer, intent(in) :: t_step real(kind(0d0)), intent(inout) :: time_avg - integer :: i, j, k, l !< Generic loop iterator real(kind(0d0)) :: start, finish call cpu_time(start) @@ -932,8 +923,6 @@ contains type(int_bounds_info) :: ix, iy, iz type(vector_field) :: gm_alpha_qp - integer :: i, j, k, l, q !< Generic loop iterator - ix%beg = 0; iy%beg = 0; iz%beg = 0 ix%end = m; iy%end = n; iz%end = p call s_convert_conservative_to_primitive_variables( & @@ -963,7 +952,7 @@ contains type(vector_field) :: gm_alpha_qp real(kind(0d0)) :: dt_local type(int_bounds_info) :: ix, iy, iz - integer :: i, j, k, l, q !< Generic loop iterators + integer :: j, k, l !< Generic loop iterators ix%beg = 0; iy%beg = 0; iz%beg = 0 ix%end = m; iy%end = n; iz%end = p