diff --git a/.cursor/rules/mfc-agent-rules.mdc b/.cursor/rules/mfc-agent-rules.mdc index a63276f493..ea67445ac4 100644 --- a/.cursor/rules/mfc-agent-rules.mdc +++ b/.cursor/rules/mfc-agent-rules.mdc @@ -1,7 +1,7 @@ ----- --description: Full MFC project rules – consolidated for Agent Mode --alwaysApply: true ----- +--- +description: Full MFC project rules – consolidated for Agent Mode +alwaysApply: true +--- # 0 Purpose & Scope Consolidated guidance for the MFC exascale, many-physics solver. @@ -19,7 +19,7 @@ Written primarily for Fortran/Fypp; the GPU and style sections matter only when - Assume free-form Fortran 2008+, `implicit none`, explicit `intent`, and modern intrinsics. - Prefer `module … contains … subroutine foo()`; avoid `COMMON` blocks and file-level `include` files. - **Read the full codebase and docs *before* changing code.** - - Docs: and the respository root `README.md`. + - Docs: and the repository root `README.md`. ### Incremental-change workflow @@ -62,27 +62,7 @@ Written primarily for Fortran/Fypp; the GPU and style sections matter only when --- -# 3 FYPP Macros for GPU acceleration Pogramming Guidelines (for GPU kernels) - -Do not directly use OpenACC or OpenMP directives directly. -Instead, use the FYPP macros contained in src/common/include/parallel_macros.fpp - -Wrap tight loops with - -```fortran -$:GPU_PARALLEL_FOR(private='[...]', copy='[...]') -``` -* Add `collapse=n` to merge nested loops when safe. -* Declare loop-local variables with `private='[...]'`. -* Allocate large arrays with `managed` or move them into a persistent - `$:GPU_ENTER_DATA(...)` region at start-up. -* **Do not** place `stop` / `error stop` inside device code. -* Must compile with Cray `ftn` and NVIDIA `nvfortran` for GPU offloading; also build CPU-only with - GNU `gfortran` and Intel `ifx`/`ifort`. - ---- - -# 4 File & Module Structure +# 3 File & Module Structure - **File Naming**: - `.fpp` files: Fypp preprocessed files that get translated to `.f90` @@ -99,25 +79,44 @@ $:GPU_PARALLEL_FOR(private='[...]', copy='[...]') - `contains` section - Implementation of subroutines and functions -# 5 Fypp Macros and GPU Acceleration +--- + +# 4 Fypp Macros -## Use of Fypp - **Fypp Directives**: - Start with `#:` (e.g., `#:include`, `#:def`, `#:enddef`) - Macros defined in `include/*.fpp` files - Used for code generation, conditional compilation, and GPU offloading -## Some examples +--- -Documentation on how to use the Fypp macros for GPU offloading is available at https://mflowcode.github.io/documentation/md_gpuParallelization.html +# 5 FYPP Macros for GPU Acceleration Programming Guidelines (for GPU kernels) -Some examples include: -- `$:GPU_ROUTINE(parallelism='[seq]')` - Marks GPU-callable routines -- `$:GPU_PARALLEL_LOOP(collapse=N)` - Parallelizes loops -- `$:GPU_LOOP(parallelism='[seq]')` - Marks sequential loops -- `$:GPU_UPDATE(device='[var1,var2]')` - Updates device data -- `$:GPU_ENTER_DATA(copyin='[var]')` - Copies data to device -- `$:GPU_EXIT_DATA(delete='[var]')` - Removes data from device +- Do not use OpenACC or OpenMP directives directly. +- Instead, use the FYPP macros contained in `src/common/include/parallel_macros.fpp` +- Documentation on how to use the Fypp macros for GPU offloading is available at https://mflowcode.github.io/documentation/md_gpuParallelization.html + +Wrap tight loops with +```fortran +$:GPU_PARALLEL_FOR(private='[...]', copy='[...]') +``` +* Add `collapse=n` to merge nested loops when safe. +* Declare loop-local variables with `private='[...]'`. +* Allocate large arrays with `managed` or move them into a persistent + `$:GPU_ENTER_DATA(...)` region at start-up. +* **Do not** place `stop` / `error stop` inside device code. +* Must compile with Cray `ftn` or NVIDIA `nvfortran` for GPU offloading; also build CPU-only with + GNU `gfortran` and Intel `ifx`/`ifort`. + +- Example GPU macros include the below, among others: + - `$:GPU_ROUTINE(parallelism='[seq]')` - Marks GPU-callable routines + - `$:GPU_PARALLEL_LOOP(collapse=N)` - Parallelizes loops + - `$:GPU_LOOP(parallelism='[seq]')` - Marks sequential loops + - `$:GPU_UPDATE(device='[var1,var2]')` - Updates device data + - `$:GPU_ENTER_DATA(copyin='[var]')` - Copies data to device + - `$:GPU_EXIT_DATA(delete='[var]')` - Removes data from device + +--- # 6 Documentation Style @@ -136,7 +135,7 @@ which conforms to the Doxygen Fortran format. - Example: `@:ASSERT(predicate, message)` - **Error Reporting**: - - Use `s_mpi_abort()` for error termination, not `stop` + - Use `s_mpi_abort(error_message)` for error termination, not `stop` - No `stop` / `error stop` inside device code # 8 Memory Management diff --git a/.fortls.json b/.fortls.json deleted file mode 100644 index a4fd1e84ad..0000000000 --- a/.fortls.json +++ /dev/null @@ -1,94 +0,0 @@ -{ - "source_dirs": [ - "src/", - "src/common/", - "src/simulation/", - "src/pre_process/", - "src/post_process/" - ], - "excl_paths": [ - "benchmarks/", - "examples/", - "tests/", - "misc/", - "src/pre_process/include/2dHardcodedIC.fpp", - "src/pre_process/include/3dHardcodedIC.fpp", - "src/pre_process/include/ExtrusionHardcodedIC.fpp" - ], - "include_dirs": [ - "src/common/include/", - "src/simulation/include/", - "src/pre_process/include/", - "src/post_process/include/" - ], - "pp_suffixes": [".fpp"], - "pp_defs": { - "MFC": 1, - "MFC_DOUBLE_PRECISION": 1 - }, - "lowercase_intrinsics": true, - "debug_log": false, - "disable_diagnostics": true, - "use_signature_help": true, - "variable_hover": true, - "hover_signature": true, - "enable_code_actions": true, - "mod_dirs": [ - "build/pre_process/", - "build/simulation/", - "build/post_process/", - "build/common/" - ], - "ext_mod_dirs": [ - "/usr/include/", - "/usr/local/include/", - "/opt/homebrew/include/" - ], - "implicit_external_mods": [ - "mpi", - "m_thermochem", - "hipfort", - "hipfort_check", - "hipfort_hipfft", - "cutensorex", - "silo_f9x", - "m_model" - ], - "disable_diagnostics_for_external_modules": true, - "max_line_length": 132, - "symbol_skip_mem": [ - "mpi_*" - ], - "disable_var_diagnostics": true, - "disable_fypp": false, - "fypp_strict": false, - "error_suppression_list": [ - "include-not-found", - "mod-not-found", - "var-masking", - "declared-twice", - "no-matching-declaration", - "invalid-parent", - "parsing-error", - "fypp-error", - "preprocessor-error", - "syntax-error", - "semantic-error", - "type-error", - "undefined-variable", - "line-too-long" - ], - "incremental_sync": false, - "debug_parser": false, - "skip_parse_errors": true, - "disable_parser": [ - "src/post_process/m_data_output.fpp", - "src/pre_process/include/ExtrusionHardcodedIC.fpp", - "src/pre_process/m_checker.fpp", - "src/pre_process/include/2dHardcodedIC.fpp", - "src/pre_process/include/3dHardcodedIC.fpp", - "src/simulation/m_qbmm.fpp", - "src/common/m_variables_conversion.fpp", - "src/simulation/m_global_parameters.fpp" - ] -} \ No newline at end of file diff --git a/.fortlsrc b/.fortlsrc index 0caf9c793a..b6b0a10bc5 100644 --- a/.fortlsrc +++ b/.fortlsrc @@ -26,7 +26,9 @@ "pp_suffixes": [".fpp"], "pp_defs": { "MFC": 1, - "MFC_DOUBLE_PRECISION": 1 + "MFC_SINGLE_PRECISION": 1, + "MFC_OPENACC": 1, + "MFC_MPI": 1 }, "lowercase_intrinsics": true, "debug_log": false, @@ -60,26 +62,11 @@ "disable_diagnostics_for_external_modules": true, "max_line_length": -1, "max_comment_line_length": -1, - "symbol_skip_mem": [ - "mpi_*" - ], "disable_var_diagnostics": false, "disable_fypp": false, "fypp_strict": false, - "error_suppression_list": [ - "include-not-found", - "mod-not-found", - "module-not-found", - "declared-twice", - "no-matching-declaration", - "invalid-parent", - "parsing-error", - "fypp-error", - "preprocessor-error", - "implicit-type" - ], "incremental_sync": false, - "debug_parser": false, + "debug_parser": true, "skip_parse_errors": true, "disable_parser": [ "src/post_process/m_data_output.fpp", diff --git a/src/common/m_boundary_common.fpp b/src/common/m_boundary_common.fpp index 66991dec10..3dc4da2c2e 100644 --- a/src/common/m_boundary_common.fpp +++ b/src/common/m_boundary_common.fpp @@ -79,17 +79,17 @@ contains !> The purpose of this procedure is to populate the buffers !! of the primitive variables, depending on the selected !! boundary conditions. - impure subroutine s_populate_variables_buffers(bc_type, q_prim_vf, pb, mv) + impure subroutine s_populate_variables_buffers(bc_type, q_prim_vf, pb_in, mv_in) type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf - real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv + real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb_in, mv_in type(integer_field), dimension(1:num_dims, -1:1), intent(in) :: bc_type integer :: k, l ! Population of Buffers in x-direction if (bc_x%beg >= 0) then - call s_mpi_sendrecv_variables_buffers(q_prim_vf, 1, -1, sys_size, pb, mv) + call s_mpi_sendrecv_variables_buffers(q_prim_vf, 1, -1, sys_size, pb_in, mv_in) else $:GPU_PARALLEL_LOOP(collapse=2) do l = 0, p @@ -98,9 +98,9 @@ contains case (BC_CHAR_SUP_OUTFLOW:BC_GHOST_EXTRAP) call s_ghost_cell_extrapolation(q_prim_vf, 1, -1, k, l) case (BC_REFLECTIVE) - call s_symmetry(q_prim_vf, 1, -1, k, l, pb, mv) + call s_symmetry(q_prim_vf, 1, -1, k, l, pb_in, mv_in) case (BC_PERIODIC) - call s_periodic(q_prim_vf, 1, -1, k, l, pb, mv) + call s_periodic(q_prim_vf, 1, -1, k, l, pb_in, mv_in) case (BC_SLIP_WALL) call s_slip_wall(q_prim_vf, 1, -1, k, l) case (BC_NO_SLIP_WALL) @@ -111,14 +111,14 @@ contains if (qbmm .and. (.not. polytropic) .and. & (bc_type(1, -1)%sf(0, k, l) <= BC_GHOST_EXTRAP)) then - call s_qbmm_extrapolation(1, -1, k, l, pb, mv) + call s_qbmm_extrapolation(1, -1, k, l, pb_in, mv_in) end if end do end do end if if (bc_x%end >= 0) then - call s_mpi_sendrecv_variables_buffers(q_prim_vf, 1, 1, sys_size, pb, mv) + call s_mpi_sendrecv_variables_buffers(q_prim_vf, 1, 1, sys_size, pb_in, mv_in) else $:GPU_PARALLEL_LOOP(collapse=2) do l = 0, p @@ -127,9 +127,9 @@ contains case (BC_CHAR_SUP_OUTFLOW:BC_GHOST_EXTRAP) ! Ghost-cell extrap. BC at end call s_ghost_cell_extrapolation(q_prim_vf, 1, 1, k, l) case (BC_REFLECTIVE) - call s_symmetry(q_prim_vf, 1, 1, k, l, pb, mv) + call s_symmetry(q_prim_vf, 1, 1, k, l, pb_in, mv_in) case (BC_PERIODIC) - call s_periodic(q_prim_vf, 1, 1, k, l, pb, mv) + call s_periodic(q_prim_vf, 1, 1, k, l, pb_in, mv_in) case (BC_SLIP_WALL) call s_slip_wall(q_prim_vf, 1, 1, k, l) case (BC_NO_SLIP_WALL) @@ -140,7 +140,7 @@ contains if (qbmm .and. (.not. polytropic) .and. & (bc_type(1, 1)%sf(0, k, l) <= BC_GHOST_EXTRAP)) then - call s_qbmm_extrapolation(1, 1, k, l, pb, mv) + call s_qbmm_extrapolation(1, 1, k, l, pb_in, mv_in) end if end do end do @@ -151,7 +151,7 @@ contains if (n == 0) return if (bc_y%beg >= 0) then - call s_mpi_sendrecv_variables_buffers(q_prim_vf, 2, -1, sys_size, pb, mv) + call s_mpi_sendrecv_variables_buffers(q_prim_vf, 2, -1, sys_size, pb_in, mv_in) else $:GPU_PARALLEL_LOOP(collapse=2) do l = 0, p @@ -160,11 +160,11 @@ contains case (BC_CHAR_SUP_OUTFLOW:BC_GHOST_EXTRAP) call s_ghost_cell_extrapolation(q_prim_vf, 2, -1, k, l) case (BC_AXIS) - call s_axis(q_prim_vf, pb, mv, k, l) + call s_axis(q_prim_vf, pb_in, mv_in, k, l) case (BC_REFLECTIVE) - call s_symmetry(q_prim_vf, 2, -1, k, l, pb, mv) + call s_symmetry(q_prim_vf, 2, -1, k, l, pb_in, mv_in) case (BC_PERIODIC) - call s_periodic(q_prim_vf, 2, -1, k, l, pb, mv) + call s_periodic(q_prim_vf, 2, -1, k, l, pb_in, mv_in) case (BC_SLIP_WALL) call s_slip_wall(q_prim_vf, 2, -1, k, l) case (BC_NO_SLIP_WALL) @@ -176,14 +176,14 @@ contains if (qbmm .and. (.not. polytropic) .and. & (bc_type(2, -1)%sf(k, 0, l) <= BC_GHOST_EXTRAP) .and. & (bc_type(2, -1)%sf(k, 0, l) /= BC_AXIS)) then - call s_qbmm_extrapolation(2, -1, k, l, pb, mv) + call s_qbmm_extrapolation(2, -1, k, l, pb_in, mv_in) end if end do end do end if if (bc_y%end >= 0) then - call s_mpi_sendrecv_variables_buffers(q_prim_vf, 2, 1, sys_size, pb, mv) + call s_mpi_sendrecv_variables_buffers(q_prim_vf, 2, 1, sys_size, pb_in, mv_in) else $:GPU_PARALLEL_LOOP(collapse=2) do l = 0, p @@ -192,9 +192,9 @@ contains case (BC_CHAR_SUP_OUTFLOW:BC_GHOST_EXTRAP) call s_ghost_cell_extrapolation(q_prim_vf, 2, 1, k, l) case (BC_REFLECTIVE) - call s_symmetry(q_prim_vf, 2, 1, k, l, pb, mv) + call s_symmetry(q_prim_vf, 2, 1, k, l, pb_in, mv_in) case (BC_PERIODIC) - call s_periodic(q_prim_vf, 2, 1, k, l, pb, mv) + call s_periodic(q_prim_vf, 2, 1, k, l, pb_in, mv_in) case (BC_SLIP_WALL) call s_slip_wall(q_prim_vf, 2, 1, k, l) case (BC_NO_SLIP_WALL) @@ -205,7 +205,7 @@ contains if (qbmm .and. (.not. polytropic) .and. & (bc_type(2, 1)%sf(k, 0, l) <= BC_GHOST_EXTRAP)) then - call s_qbmm_extrapolation(2, 1, k, l, pb, mv) + call s_qbmm_extrapolation(2, 1, k, l, pb_in, mv_in) end if end do end do @@ -216,7 +216,7 @@ contains if (p == 0) return if (bc_z%beg >= 0) then - call s_mpi_sendrecv_variables_buffers(q_prim_vf, 3, -1, sys_size, pb, mv) + call s_mpi_sendrecv_variables_buffers(q_prim_vf, 3, -1, sys_size, pb_in, mv_in) else $:GPU_PARALLEL_LOOP(collapse=2) do l = -buff_size, n + buff_size @@ -225,9 +225,9 @@ contains case (BC_CHAR_SUP_OUTFLOW:BC_GHOST_EXTRAP) call s_ghost_cell_extrapolation(q_prim_vf, 3, -1, k, l) case (BC_REFLECTIVE) - call s_symmetry(q_prim_vf, 3, -1, k, l, pb, mv) + call s_symmetry(q_prim_vf, 3, -1, k, l, pb_in, mv_in) case (BC_PERIODIC) - call s_periodic(q_prim_vf, 3, -1, k, l, pb, mv) + call s_periodic(q_prim_vf, 3, -1, k, l, pb_in, mv_in) case (BC_SLIP_WALL) call s_slip_wall(q_prim_vf, 3, -1, k, l) case (BC_NO_SLIP_WALL) @@ -238,14 +238,14 @@ contains if (qbmm .and. (.not. polytropic) .and. & (bc_type(3, -1)%sf(k, l, 0) <= BC_GHOST_EXTRAP)) then - call s_qbmm_extrapolation(3, -1, k, l, pb, mv) + call s_qbmm_extrapolation(3, -1, k, l, pb_in, mv_in) end if end do end do end if if (bc_z%end >= 0) then - call s_mpi_sendrecv_variables_buffers(q_prim_vf, 3, 1, sys_size, pb, mv) + call s_mpi_sendrecv_variables_buffers(q_prim_vf, 3, 1, sys_size, pb_in, mv_in) else $:GPU_PARALLEL_LOOP(collapse=2) do l = -buff_size, n + buff_size @@ -254,9 +254,9 @@ contains case (BC_CHAR_SUP_OUTFLOW:BC_GHOST_EXTRAP) call s_ghost_cell_extrapolation(q_prim_vf, 3, 1, k, l) case (BC_REFLECTIVE) - call s_symmetry(q_prim_vf, 3, 1, k, l, pb, mv) + call s_symmetry(q_prim_vf, 3, 1, k, l, pb_in, mv_in) case (BC_PERIODIC) - call s_periodic(q_prim_vf, 3, 1, k, l, pb, mv) + call s_periodic(q_prim_vf, 3, 1, k, l, pb_in, mv_in) case (BC_SlIP_WALL) call s_slip_wall(q_prim_vf, 3, 1, k, l) case (BC_NO_SLIP_WALL) @@ -267,7 +267,7 @@ contains if (qbmm .and. (.not. polytropic) .and. & (bc_type(3, 1)%sf(k, l, 0) <= BC_GHOST_EXTRAP)) then - call s_qbmm_extrapolation(3, 1, k, l, pb, mv) + call s_qbmm_extrapolation(3, 1, k, l, pb_in, mv_in) end if end do end do @@ -337,10 +337,10 @@ contains end subroutine s_ghost_cell_extrapolation - pure subroutine s_symmetry(q_prim_vf, bc_dir, bc_loc, k, l, pb, mv) + pure subroutine s_symmetry(q_prim_vf, bc_dir, bc_loc, k, l, pb_in, mv_in) $:GPU_ROUTINE(parallelism='[seq]') type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf - real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv + real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb_in, mv_in integer, intent(in) :: bc_dir, bc_loc integer, intent(in) :: k, l @@ -380,10 +380,10 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(-j, k, l, q, i) = & - pb(j - 1, k, l, q, i) - mv(-j, k, l, q, i) = & - mv(j - 1, k, l, q, i) + pb_in(-j, k, l, q, i) = & + pb_in(j - 1, k, l, q, i) + mv_in(-j, k, l, q, i) = & + mv_in(j - 1, k, l, q, i) end do end do end do @@ -420,10 +420,10 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(m + j, k, l, q, i) = & - pb(m - (j - 1), k, l, q, i) - mv(m + j, k, l, q, i) = & - mv(m - (j - 1), k, l, q, i) + pb_in(m + j, k, l, q, i) = & + pb_in(m - (j - 1), k, l, q, i) + mv_in(m + j, k, l, q, i) = & + mv_in(m - (j - 1), k, l, q, i) end do end do end do @@ -462,10 +462,10 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(k, -j, l, q, i) = & - pb(k, j - 1, l, q, i) - mv(k, -j, l, q, i) = & - mv(k, j - 1, l, q, i) + pb_in(k, -j, l, q, i) = & + pb_in(k, j - 1, l, q, i) + mv_in(k, -j, l, q, i) = & + mv_in(k, j - 1, l, q, i) end do end do end do @@ -502,10 +502,10 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(k, n + j, l, q, i) = & - pb(k, n - (j - 1), l, q, i) - mv(k, n + j, l, q, i) = & - mv(k, n - (j - 1), l, q, i) + pb_in(k, n + j, l, q, i) = & + pb_in(k, n - (j - 1), l, q, i) + mv_in(k, n + j, l, q, i) = & + mv_in(k, n - (j - 1), l, q, i) end do end do end do @@ -544,10 +544,10 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(k, l, -j, q, i) = & - pb(k, l, j - 1, q, i) - mv(k, l, -j, q, i) = & - mv(k, l, j - 1, q, i) + pb_in(k, l, -j, q, i) = & + pb_in(k, l, j - 1, q, i) + mv_in(k, l, -j, q, i) = & + mv_in(k, l, j - 1, q, i) end do end do end do @@ -584,10 +584,10 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(k, l, p + j, q, i) = & - pb(k, l, p - (j - 1), q, i) - mv(k, l, p + j, q, i) = & - mv(k, l, p - (j - 1), q, i) + pb_in(k, l, p + j, q, i) = & + pb_in(k, l, p - (j - 1), q, i) + mv_in(k, l, p + j, q, i) = & + mv_in(k, l, p - (j - 1), q, i) end do end do end do @@ -597,10 +597,10 @@ contains end subroutine s_symmetry - pure subroutine s_periodic(q_prim_vf, bc_dir, bc_loc, k, l, pb, mv) + pure subroutine s_periodic(q_prim_vf, bc_dir, bc_loc, k, l, pb_in, mv_in) $:GPU_ROUTINE(parallelism='[seq]') type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf - real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv + real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb_in, mv_in integer, intent(in) :: bc_dir, bc_loc integer, intent(in) :: k, l @@ -619,10 +619,10 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(-j, k, l, q, i) = & - pb(m - (j - 1), k, l, q, i) - mv(-j, k, l, q, i) = & - mv(m - (j - 1), k, l, q, i) + pb_in(-j, k, l, q, i) = & + pb_in(m - (j - 1), k, l, q, i) + mv_in(-j, k, l, q, i) = & + mv_in(m - (j - 1), k, l, q, i) end do end do end do @@ -639,10 +639,10 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(m + j, k, l, q, i) = & - pb(j - 1, k, l, q, i) - mv(m + j, k, l, q, i) = & - mv(j - 1, k, l, q, i) + pb_in(m + j, k, l, q, i) = & + pb_in(j - 1, k, l, q, i) + mv_in(m + j, k, l, q, i) = & + mv_in(j - 1, k, l, q, i) end do end do end do @@ -661,10 +661,10 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(k, -j, l, q, i) = & - pb(k, n - (j - 1), l, q, i) - mv(k, -j, l, q, i) = & - mv(k, n - (j - 1), l, q, i) + pb_in(k, -j, l, q, i) = & + pb_in(k, n - (j - 1), l, q, i) + mv_in(k, -j, l, q, i) = & + mv_in(k, n - (j - 1), l, q, i) end do end do end do @@ -681,10 +681,10 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(k, n + j, l, q, i) = & - pb(k, (j - 1), l, q, i) - mv(k, n + j, l, q, i) = & - mv(k, (j - 1), l, q, i) + pb_in(k, n + j, l, q, i) = & + pb_in(k, (j - 1), l, q, i) + mv_in(k, n + j, l, q, i) = & + mv_in(k, (j - 1), l, q, i) end do end do end do @@ -703,10 +703,10 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(k, l, -j, q, i) = & - pb(k, l, p - (j - 1), q, i) - mv(k, l, -j, q, i) = & - mv(k, l, p - (j - 1), q, i) + pb_in(k, l, -j, q, i) = & + pb_in(k, l, p - (j - 1), q, i) + mv_in(k, l, -j, q, i) = & + mv_in(k, l, p - (j - 1), q, i) end do end do end do @@ -723,10 +723,10 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(k, l, p + j, q, i) = & - pb(k, l, j - 1, q, i) - mv(k, l, p + j, q, i) = & - mv(k, l, j - 1, q, i) + pb_in(k, l, p + j, q, i) = & + pb_in(k, l, j - 1, q, i) + mv_in(k, l, p + j, q, i) = & + mv_in(k, l, j - 1, q, i) end do end do end do @@ -736,10 +736,10 @@ contains end subroutine s_periodic - pure subroutine s_axis(q_prim_vf, pb, mv, k, l) + pure subroutine s_axis(q_prim_vf, pb_in, mv_in, k, l) $:GPU_ROUTINE(parallelism='[seq]') type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf - real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv + real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb_in, mv_in integer, intent(in) :: k, l integer :: j, q, i @@ -784,10 +784,10 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(k, -j, l, q, i) = & - pb(k, j - 1, l - ((p + 1)/2), q, i) - mv(k, -j, l, q, i) = & - mv(k, j - 1, l - ((p + 1)/2), q, i) + pb_in(k, -j, l, q, i) = & + pb_in(k, j - 1, l - ((p + 1)/2), q, i) + mv_in(k, -j, l, q, i) = & + mv_in(k, j - 1, l - ((p + 1)/2), q, i) end do end do end do @@ -1079,9 +1079,9 @@ contains end subroutine s_dirichlet - pure subroutine s_qbmm_extrapolation(bc_dir, bc_loc, k, l, pb, mv) + pure subroutine s_qbmm_extrapolation(bc_dir, bc_loc, k, l, pb_in, mv_in) $:GPU_ROUTINE(parallelism='[seq]') - real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv + real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb_in, mv_in integer, intent(in) :: bc_dir, bc_loc integer, intent(in) :: k, l @@ -1092,8 +1092,8 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(-j, k, l, q, i) = pb(0, k, l, q, i) - mv(-j, k, l, q, i) = mv(0, k, l, q, i) + pb_in(-j, k, l, q, i) = pb_in(0, k, l, q, i) + mv_in(-j, k, l, q, i) = mv_in(0, k, l, q, i) end do end do end do @@ -1101,8 +1101,8 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(m + j, k, l, q, i) = pb(m, k, l, q, i) - mv(m + j, k, l, q, i) = mv(m, k, l, q, i) + pb_in(m + j, k, l, q, i) = pb_in(m, k, l, q, i) + mv_in(m + j, k, l, q, i) = mv_in(m, k, l, q, i) end do end do end do @@ -1112,8 +1112,8 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(k, -j, l, q, i) = pb(k, 0, l, q, i) - mv(k, -j, l, q, i) = mv(k, 0, l, q, i) + pb_in(k, -j, l, q, i) = pb_in(k, 0, l, q, i) + mv_in(k, -j, l, q, i) = mv_in(k, 0, l, q, i) end do end do end do @@ -1121,8 +1121,8 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(k, n + j, l, q, i) = pb(k, n, l, q, i) - mv(k, n + j, l, q, i) = mv(k, n, l, q, i) + pb_in(k, n + j, l, q, i) = pb_in(k, n, l, q, i) + mv_in(k, n + j, l, q, i) = mv_in(k, n, l, q, i) end do end do end do @@ -1132,8 +1132,8 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(k, l, -j, q, i) = pb(k, l, 0, q, i) - mv(k, l, -j, q, i) = mv(k, l, 0, q, i) + pb_in(k, l, -j, q, i) = pb_in(k, l, 0, q, i) + mv_in(k, l, -j, q, i) = mv_in(k, l, 0, q, i) end do end do end do @@ -1141,8 +1141,8 @@ contains do i = 1, nb do q = 1, nnode do j = 1, buff_size - pb(k, l, p + j, q, i) = pb(k, l, p, q, i) - mv(k, l, p + j, q, i) = mv(k, l, p, q, i) + pb_in(k, l, p + j, q, i) = pb_in(k, l, p, q, i) + mv_in(k, l, p + j, q, i) = mv_in(k, l, p, q, i) end do end do end do @@ -1657,11 +1657,11 @@ contains #endif end subroutine s_create_mpi_types - subroutine s_write_serial_boundary_condition_files(q_prim_vf, bc_type, step_dirpath, old_grid) + subroutine s_write_serial_boundary_condition_files(q_prim_vf, bc_type, step_dirpath, old_grid_in) type(scalar_field), dimension(sys_size) :: q_prim_vf type(integer_field), dimension(1:num_dims, -1:1) :: bc_type - logical :: old_grid + logical :: old_grid_in character(LEN=*), intent(in) :: step_dirpath @@ -1670,7 +1670,7 @@ contains character(len=10) :: status - if (old_grid) then + if (old_grid_in) then status = 'old' else status = 'new' diff --git a/src/common/m_finite_differences.fpp b/src/common/m_finite_differences.fpp index 1857a31cd8..2430374f4f 100644 --- a/src/common/m_finite_differences.fpp +++ b/src/common/m_finite_differences.fpp @@ -69,17 +69,17 @@ contains !! @param q Number of cells in the s-coordinate direction !! @param s_cc Locations of the cell-centers in the s-coordinate direction !! @param fd_coeff_s Finite-diff. coefficients in the s-coordinate direction - pure subroutine s_compute_finite_difference_coefficients(q, s_cc, fd_coeff_s, buff_size, & + pure subroutine s_compute_finite_difference_coefficients(q, s_cc, fd_coeff_s, local_buff_size, & fd_number_in, fd_order_in, offset_s) integer :: lB, lE !< loop bounds integer, intent(IN) :: q - integer, intent(IN) :: buff_size, fd_number_in, fd_order_in + integer, intent(IN) :: local_buff_size, fd_number_in, fd_order_in type(int_bounds_info), optional, intent(IN) :: offset_s real(wp), allocatable, dimension(:, :), intent(INOUT) :: fd_coeff_s real(wp), & - dimension(-buff_size:q + buff_size), & + dimension(-local_buff_size:q + local_buff_size), & intent(IN) :: s_cc integer :: i !< Generic loop iterator diff --git a/src/common/m_helper.fpp b/src/common/m_helper.fpp index a067853c66..f432ca12fa 100644 --- a/src/common/m_helper.fpp +++ b/src/common/m_helper.fpp @@ -79,11 +79,11 @@ contains real(wp), optional, intent(in) :: div integer :: i, j - integer :: m, n + integer :: local_m, local_n real(wp) :: c - m = size(A, 1) - n = size(A, 2) + local_m = size(A, 1) + local_n = size(A, 2) if (present(div)) then c = div @@ -91,10 +91,10 @@ contains c = 1._wp end if - print *, m, n + print *, local_m, local_n - do i = 1, m - do j = 1, n + do i = 1, local_m + do j = 1, local_n write (*, fmt="(F12.4)", advance="no") A(i, j)/c end do write (*, fmt="(A1)") " " @@ -218,7 +218,7 @@ contains end subroutine s_initialize_nonpoly !> Computes the transfer coefficient for the non-polytropic bubble compression process - !! @param omega natural frqeuencies + !! @param omega natural frequencies !! @param peclet Peclet number !! @param Re_trans Real part of the transport coefficients !! @param Im_trans Imaginary part of the transport coefficients @@ -251,10 +251,10 @@ contains end subroutine s_int_to_str !> Computes the Simpson weights for quadrature - subroutine s_simpson(weight, R0) + subroutine s_simpson(local_weight, local_R0) - real(wp), dimension(:), intent(inout) :: weight - real(wp), dimension(:), intent(inout) :: R0 + real(wp), dimension(:), intent(inout) :: local_weight + real(wp), dimension(:), intent(inout) :: local_R0 integer :: ir real(wp) :: R0mn, R0mx, dphi, tmp, sd @@ -268,7 +268,7 @@ contains do ir = 1, nb phi(ir) = log(R0mn) & + (ir - 1._wp)*log(R0mx/R0mn)/(nb - 1._wp) - R0(ir) = exp(phi(ir)) + local_R0(ir) = exp(phi(ir)) end do dphi = phi(2) - phi(1) @@ -277,15 +277,15 @@ contains ! Gaussian tmp = exp(-0.5_wp*(phi(ir)/sd)**2)/sqrt(2._wp*pi)/sd if (mod(ir, 2) == 0) then - weight(ir) = tmp*4._wp*dphi/3._wp + local_weight(ir) = tmp*4._wp*dphi/3._wp else - weight(ir) = tmp*2._wp*dphi/3._wp + local_weight(ir) = tmp*2._wp*dphi/3._wp end if end do tmp = exp(-0.5_wp*(phi(1)/sd)**2)/sqrt(2._wp*pi)/sd - weight(1) = tmp*dphi/3._wp + local_weight(1) = tmp*dphi/3._wp tmp = exp(-0.5_wp*(phi(nb)/sd)**2)/sqrt(2._wp*pi)/sd - weight(nb) = tmp*dphi/3._wp + local_weight(nb) = tmp*dphi/3._wp end subroutine s_simpson !> This procedure computes the cross product of two vectors. @@ -318,40 +318,40 @@ contains !> This procedure creates a transformation matrix. !! @param p Parameters for the transformation. !! @return Transformation matrix. - pure function f_create_transform_matrix(p, center) result(out_matrix) + pure function f_create_transform_matrix(param, center) result(out_matrix) - type(ic_model_parameters), intent(in) :: p + type(ic_model_parameters), intent(in) :: param real(wp), dimension(1:3), optional, intent(in) :: center real(wp), dimension(1:4, 1:4) :: sc, rz, rx, ry, tr, t_back, t_to_origin, out_matrix sc = transpose(reshape([ & - p%scale(1), 0._wp, 0._wp, 0._wp, & - 0._wp, p%scale(2), 0._wp, 0._wp, & - 0._wp, 0._wp, p%scale(3), 0._wp, & + param%scale(1), 0._wp, 0._wp, 0._wp, & + 0._wp, param%scale(2), 0._wp, 0._wp, & + 0._wp, 0._wp, param%scale(3), 0._wp, & 0._wp, 0._wp, 0._wp, 1._wp], shape(sc))) rz = transpose(reshape([ & - cos(p%rotate(3)), -sin(p%rotate(3)), 0._wp, 0._wp, & - sin(p%rotate(3)), cos(p%rotate(3)), 0._wp, 0._wp, & + cos(param%rotate(3)), -sin(param%rotate(3)), 0._wp, 0._wp, & + sin(param%rotate(3)), cos(param%rotate(3)), 0._wp, 0._wp, & 0._wp, 0._wp, 1._wp, 0._wp, & 0._wp, 0._wp, 0._wp, 1._wp], shape(rz))) rx = transpose(reshape([ & 1._wp, 0._wp, 0._wp, 0._wp, & - 0._wp, cos(p%rotate(1)), -sin(p%rotate(1)), 0._wp, & - 0._wp, sin(p%rotate(1)), cos(p%rotate(1)), 0._wp, & + 0._wp, cos(param%rotate(1)), -sin(param%rotate(1)), 0._wp, & + 0._wp, sin(param%rotate(1)), cos(param%rotate(1)), 0._wp, & 0._wp, 0._wp, 0._wp, 1._wp], shape(rx))) ry = transpose(reshape([ & - cos(p%rotate(2)), 0._wp, sin(p%rotate(2)), 0._wp, & + cos(param%rotate(2)), 0._wp, sin(param%rotate(2)), 0._wp, & 0._wp, 1._wp, 0._wp, 0._wp, & - -sin(p%rotate(2)), 0._wp, cos(p%rotate(2)), 0._wp, & + -sin(param%rotate(2)), 0._wp, cos(param%rotate(2)), 0._wp, & 0._wp, 0._wp, 0._wp, 1._wp], shape(ry))) tr = transpose(reshape([ & - 1._wp, 0._wp, 0._wp, p%translate(1), & - 0._wp, 1._wp, 0._wp, p%translate(2), & - 0._wp, 0._wp, 1._wp, p%translate(3), & + 1._wp, 0._wp, 0._wp, param%translate(1), & + 0._wp, 1._wp, 0._wp, param%translate(2), & + 0._wp, 0._wp, 1._wp, param%translate(3), & 0._wp, 0._wp, 0._wp, 1._wp], shape(tr))) if (present(center)) then @@ -484,18 +484,18 @@ contains !! @param x is the input value !! @param l is the degree !! @return P is the unassociated legendre polynomial evaluated at x - pure recursive function unassociated_legendre(x, l) result(P) + pure recursive function unassociated_legendre(x, l) result(result_P) integer, intent(in) :: l real(wp), intent(in) :: x - real(wp) :: P + real(wp) :: result_P if (l == 0) then - P = 1._wp + result_P = 1._wp else if (l == 1) then - P = x + result_P = x else - P = ((2*l - 1)*x*unassociated_legendre(x, l - 1) - (l - 1)*unassociated_legendre(x, l - 2))/l + result_P = ((2*l - 1)*x*unassociated_legendre(x, l - 1) - (l - 1)*unassociated_legendre(x, l - 2))/l end if end function unassociated_legendre @@ -504,20 +504,20 @@ contains !! @param x is the x coordinate !! @param phi is the phi coordinate !! @param l is the degree - !! @param m is the order + !! @param m_order is the order !! @return Y is the spherical harmonic function evaluated at x and phi - pure recursive function spherical_harmonic_func(x, phi, l, m) result(Y) + pure recursive function spherical_harmonic_func(x, phi, l, m_order) result(Y) - integer, intent(in) :: l, m + integer, intent(in) :: l, m_order real(wp), intent(in) :: x, phi - real(wp) :: Y, prefactor, pi - - pi = acos(-1._wp) - prefactor = sqrt((2*l + 1)/(4*pi)*factorial(l - m)/factorial(l + m)); - if (m == 0) then - Y = prefactor*associated_legendre(x, l, m); - elseif (m > 0) then - Y = (-1._wp)**m*sqrt(2._wp)*prefactor*associated_legendre(x, l, m)*cos(m*phi); + real(wp) :: Y, prefactor, local_pi + + local_pi = acos(-1._wp) + prefactor = sqrt((2*l + 1)/(4*local_pi)*factorial(l - m_order)/factorial(l + m_order)); + if (m_order == 0) then + Y = prefactor*associated_legendre(x, l, m_order); + elseif (m_order > 0) then + Y = (-1._wp)**m_order*sqrt(2._wp)*prefactor*associated_legendre(x, l, m_order)*cos(m_order*phi); end if end function spherical_harmonic_func @@ -526,56 +526,56 @@ contains !! at x with inputs l and m !! @param x is the input value !! @param l is the degree - !! @param m is the order + !! @param m_order is the order !! @return P is the associated legendre polynomial evaluated at x - pure recursive function associated_legendre(x, l, m) result(P) + pure recursive function associated_legendre(x, l, m_order) result(result_P) - integer, intent(in) :: l, m + integer, intent(in) :: l, m_order real(wp), intent(in) :: x - real(wp) :: P - - if (m <= 0 .and. l <= 0) then - P = 1; - elseif (l == 1 .and. m <= 0) then - P = x; - elseif (l == 1 .and. m == 1) then - P = -(1 - x**2)**(1._wp/2._wp); - elseif (m == l) then - P = (-1)**l*double_factorial(2*l - 1)*(1 - x**2)**(l/2); - elseif (m == l - 1) then - P = x*(2*l - 1)*associated_legendre(x, l - 1, l - 1); + real(wp) :: result_P + + if (m_order <= 0 .and. l <= 0) then + result_P = 1; + elseif (l == 1 .and. m_order <= 0) then + result_P = x; + elseif (l == 1 .and. m_order == 1) then + result_P = -(1 - x**2)**(1._wp/2._wp); + elseif (m_order == l) then + result_P = (-1)**l*double_factorial(2*l - 1)*(1 - x**2)**(l/2); + elseif (m_order == l - 1) then + result_P = x*(2*l - 1)*associated_legendre(x, l - 1, l - 1); else - P = ((2*l - 1)*x*associated_legendre(x, l - 1, m) - (l + m - 1)*associated_legendre(x, l - 2, m))/(l - m); + result_P = ((2*l - 1)*x*associated_legendre(x, l - 1, m_order) - (l + m_order - 1)*associated_legendre(x, l - 2, m_order))/(l - m_order); end if end function associated_legendre !> This function calculates the double factorial value of an integer - !! @param n is the input integer + !! @param n_in is the input integer !! @return R is the double factorial value of n - pure elemental function double_factorial(n) result(R) + pure elemental function double_factorial(n_in) result(R_result) - integer, intent(in) :: n + integer, intent(in) :: n_in integer, parameter :: int64_kind = selected_int_kind(18) ! 18 bytes for 64-bit integer - integer(kind=int64_kind) :: R + integer(kind=int64_kind) :: R_result integer :: i - R = product((/(i, i=n, 1, -2)/)) + R_result = product((/(i, i=n_in, 1, -2)/)) end function double_factorial !> The following function calculates the factorial value of an integer - !! @param n is the input integer + !! @param n_in is the input integer !! @return R is the factorial value of n - pure elemental function factorial(n) result(R) + pure elemental function factorial(n_in) result(R_result) - integer, intent(in) :: n + integer, intent(in) :: n_in integer, parameter :: int64_kind = selected_int_kind(18) ! 18 bytes for 64-bit integer - integer(kind=int64_kind) :: R + integer(kind=int64_kind) :: R_result integer :: i - R = product((/(i, i=n, 1, -1)/)) + R_result = product((/(i, i=n_in, 1, -1)/)) end function factorial diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index 496f2a49b9..da485fa3ac 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -24,7 +24,7 @@ module m_mpi_common implicit none - integer, private :: ierr, v_size !< + integer, private :: v_size $:GPU_DECLARE(create='[v_size]') !! Generic flags used to identify and report MPI errors @@ -88,6 +88,10 @@ contains !! available for the job and the local processor rank. impure subroutine s_mpi_initialize +#ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors +#endif + #ifndef MFC_MPI ! Serial run only has 1 processor @@ -136,6 +140,7 @@ contains ! Generic loop iterator integer :: i, j + integer :: ierr !< Generic flag used to identify and report MPI errors !Altered system size for the lagrangian subgrid bubble model integer :: alt_sys @@ -284,7 +289,8 @@ contains integer, intent(in) :: root ! Rank of the root process real(wp), allocatable, intent(out) :: gathered_vector(:) ! Gathered vector on the root process - integer :: i, ierr + integer :: i + integer :: ierr !< Generic flag used to identify and report MPI errors integer, allocatable :: recounts(:), displs(:) #ifdef MFC_MPI @@ -314,6 +320,7 @@ contains real(wp), intent(inout) :: time_avg #ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors call MPI_GATHER(time_avg, 1, mpi_p, proc_time(0), 1, mpi_p, 0, MPI_COMM_WORLD, ierr) @@ -365,6 +372,7 @@ contains #ifdef MFC_SIMULATION #ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors ! Reducing local extrema of ICFL, VCFL, CCFL and Rc numbers to their ! global extrema and bookkeeping the results on the rank 0 processor @@ -408,6 +416,7 @@ contains real(wp), intent(out) :: var_glb #ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors ! Performing the reduction procedure call MPI_ALLREDUCE(var_loc, var_glb, 1, mpi_p, & @@ -430,6 +439,7 @@ contains real(wp), intent(out) :: var_glb #ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors ! Performing the reduction procedure call MPI_ALLREDUCE(var_loc, var_glb, 1, mpi_p, & @@ -452,6 +462,7 @@ contains real(wp), intent(out) :: var_glb #ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors ! Performing the reduction procedure call MPI_ALLREDUCE(var_loc, var_glb, 1, mpi_p, & @@ -472,6 +483,7 @@ contains real(wp), intent(inout) :: var_loc #ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors ! Temporary storage variable that holds the reduced minimum value real(wp) :: var_glb @@ -507,6 +519,7 @@ contains real(wp), dimension(2), intent(inout) :: var_loc #ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors real(wp), dimension(2) :: var_glb !< !! Temporary storage variable that holds the reduced maximum value @@ -533,6 +546,10 @@ contains character(len=*), intent(in), optional :: prnt integer, intent(in), optional :: code +#ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors +#endif + if (present(prnt)) then print *, prnt call flush (6) @@ -560,6 +577,7 @@ contains impure subroutine s_mpi_barrier #ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors ! Calling MPI_BARRIER call MPI_BARRIER(MPI_COMM_WORLD, ierr) @@ -572,6 +590,7 @@ contains impure subroutine s_mpi_finalize #ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors ! Finalizing the MPI environment call MPI_FINALIZE(ierr) @@ -590,10 +609,10 @@ contains mpi_dir, & pbc_loc, & nVar, & - pb, mv) + pb_in, mv_in) type(scalar_field), dimension(1:), intent(inout) :: q_comm - real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv + real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb_in, mv_in integer, intent(in) :: mpi_dir, pbc_loc, nVar integer :: i, j, k, l, r, q !< Generic loop iterators @@ -609,12 +628,13 @@ contains integer :: pack_offset, unpack_offset #ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors call nvtxStartRange("RHS-COMM-PACKBUF") qbmm_comm = .false. - if (present(pb) .and. present(mv) .and. qbmm .and. .not. polytropic) then + if (present(pb_in) .and. present(mv_in) .and. qbmm .and. .not. polytropic) then qbmm_comm = .true. v_size = nVar + 2*nb*4 buffer_counts = (/ & @@ -688,7 +708,7 @@ contains do q = 1, nb r = (i - 1) + (q - 1)*4 + v_size* & (j + buff_size*(k + (n + 1)*l)) - buff_send(r) = pb(j + pack_offset, k, l, i - nVar, q) + buff_send(r) = pb_in(j + pack_offset, k, l, i - nVar, q) end do end do end do @@ -703,7 +723,7 @@ contains do q = 1, nb r = (i - 1) + (q - 1)*4 + nb*4 + v_size* & (j + buff_size*(k + (n + 1)*l)) - buff_send(r) = mv(j + pack_offset, k, l, i - nVar, q) + buff_send(r) = mv_in(j + pack_offset, k, l, i - nVar, q) end do end do end do @@ -735,7 +755,7 @@ contains r = (i - 1) + (q - 1)*4 + v_size* & ((j + buff_size) + (m + 2*buff_size + 1)* & (k + buff_size*l)) - buff_send(r) = pb(j, k + pack_offset, l, i - nVar, q) + buff_send(r) = pb_in(j, k + pack_offset, l, i - nVar, q) end do end do end do @@ -751,7 +771,7 @@ contains r = (i - 1) + (q - 1)*4 + nb*4 + v_size* & ((j + buff_size) + (m + 2*buff_size + 1)* & (k + buff_size*l)) - buff_send(r) = mv(j, k + pack_offset, l, i - nVar, q) + buff_send(r) = mv_in(j, k + pack_offset, l, i - nVar, q) end do end do end do @@ -783,7 +803,7 @@ contains r = (i - 1) + (q - 1)*4 + v_size* & ((j + buff_size) + (m + 2*buff_size + 1)* & ((k + buff_size) + (n + 2*buff_size + 1)*l)) - buff_send(r) = pb(j, k, l + pack_offset, i - nVar, q) + buff_send(r) = pb_in(j, k, l + pack_offset, i - nVar, q) end do end do end do @@ -799,7 +819,7 @@ contains r = (i - 1) + (q - 1)*4 + nb*4 + v_size* & ((j + buff_size) + (m + 2*buff_size + 1)* & ((k + buff_size) + (n + 2*buff_size + 1)*l)) - buff_send(r) = mv(j, k, l + pack_offset, i - nVar, q) + buff_send(r) = mv_in(j, k, l + pack_offset, i - nVar, q) end do end do end do @@ -887,7 +907,7 @@ contains do q = 1, nb r = (i - 1) + (q - 1)*4 + v_size* & (j + buff_size*((k + 1) + (n + 1)*l)) - pb(j + unpack_offset, k, l, i - nVar, q) = buff_recv(r) + pb_in(j + unpack_offset, k, l, i - nVar, q) = buff_recv(r) end do end do end do @@ -902,7 +922,7 @@ contains do q = 1, nb r = (i - 1) + (q - 1)*4 + nb*4 + v_size* & (j + buff_size*((k + 1) + (n + 1)*l)) - mv(j + unpack_offset, k, l, i - nVar, q) = buff_recv(r) + mv_in(j + unpack_offset, k, l, i - nVar, q) = buff_recv(r) end do end do end do @@ -940,7 +960,7 @@ contains r = (i - 1) + (q - 1)*4 + v_size* & ((j + buff_size) + (m + 2*buff_size + 1)* & ((k + buff_size) + buff_size*l)) - pb(j, k + unpack_offset, l, i - nVar, q) = buff_recv(r) + pb_in(j, k + unpack_offset, l, i - nVar, q) = buff_recv(r) end do end do end do @@ -956,7 +976,7 @@ contains r = (i - 1) + (q - 1)*4 + nb*4 + v_size* & ((j + buff_size) + (m + 2*buff_size + 1)* & ((k + buff_size) + buff_size*l)) - mv(j, k + unpack_offset, l, i - nVar, q) = buff_recv(r) + mv_in(j, k + unpack_offset, l, i - nVar, q) = buff_recv(r) end do end do end do @@ -997,7 +1017,7 @@ contains ((j + buff_size) + (m + 2*buff_size + 1)* & ((k + buff_size) + (n + 2*buff_size + 1)* & (l + buff_size))) - pb(j, k, l + unpack_offset, i - nVar, q) = buff_recv(r) + pb_in(j, k, l + unpack_offset, i - nVar, q) = buff_recv(r) end do end do end do @@ -1014,7 +1034,7 @@ contains ((j + buff_size) + (m + 2*buff_size + 1)* & ((k + buff_size) + (n + 2*buff_size + 1)* & (l + buff_size))) - mv(j, k, l + unpack_offset, i - nVar, q) = buff_recv(r) + mv_in(j, k, l + unpack_offset, i - nVar, q) = buff_recv(r) end do end do end do @@ -1058,6 +1078,7 @@ contains integer :: recon_order !< reconstruction order integer :: i, j !< Generic loop iterators + integer :: ierr !< Generic flag used to identify and report MPI errors if (num_procs == 1 .and. parallel_io) then do i = 1, num_dims @@ -1532,6 +1553,7 @@ contains integer, intent(in) :: pbc_loc #ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors ! MPI Communication in x-direction if (mpi_dir == 1) then diff --git a/src/post_process/m_data_input.f90 b/src/post_process/m_data_input.f90 index 8409f8776f..f0b4e5da6a 100644 --- a/src/post_process/m_data_input.f90 +++ b/src/post_process/m_data_input.f90 @@ -196,22 +196,22 @@ end subroutine s_read_ib_data_files !> Helper subroutine to allocate field arrays for given dimensionality !! @param start_idx Starting index for allocation !! @param end_x, end_y, end_z End indices for each dimension - impure subroutine s_allocate_field_arrays(start_idx, end_x, end_y, end_z) + impure subroutine s_allocate_field_arrays(local_start_idx, end_x, end_y, end_z) - integer, intent(in) :: start_idx, end_x, end_y, end_z + integer, intent(in) :: local_start_idx, end_x, end_y, end_z integer :: i do i = 1, sys_size - allocate (q_cons_vf(i)%sf(start_idx:end_x, start_idx:end_y, start_idx:end_z)) - allocate (q_prim_vf(i)%sf(start_idx:end_x, start_idx:end_y, start_idx:end_z)) + allocate (q_cons_vf(i)%sf(local_start_idx:end_x, local_start_idx:end_y, local_start_idx:end_z)) + allocate (q_prim_vf(i)%sf(local_start_idx:end_x, local_start_idx:end_y, local_start_idx:end_z)) end do if (ib) then - allocate (ib_markers%sf(start_idx:end_x, start_idx:end_y, start_idx:end_z)) + allocate (ib_markers%sf(local_start_idx:end_x, local_start_idx:end_y, local_start_idx:end_z)) end if if (chemistry) then - allocate (q_T_sf%sf(start_idx:end_x, start_idx:end_y, start_idx:end_z)) + allocate (q_T_sf%sf(local_start_idx:end_x, local_start_idx:end_y, local_start_idx:end_z)) end if end subroutine s_allocate_field_arrays diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp index a6e7b12c24..4c8867225e 100644 --- a/src/post_process/m_data_output.fpp +++ b/src/post_process/m_data_output.fpp @@ -105,7 +105,7 @@ module m_data_output ! Generic error flags utilized in the handling, checking and the reporting ! of the input and output operations errors with a formatted database file - integer, private :: err, ierr + integer, private :: err contains @@ -477,6 +477,8 @@ contains ! Generic string used to store the location of a particular file character(LEN=len_trim(case_dir) + 3*name_len) :: file_loc + integer :: ierr !< Generic flag used to identify and report database errors + ! Silo-HDF5 Database Format if (format == 1) then @@ -650,6 +652,8 @@ contains ! Generic loop iterator integer :: i + integer :: ierr !< Generic flag used to identify and report database errors + ! Silo-HDF5 Database Format if (format == 1 .and. n > 0) then @@ -860,6 +864,8 @@ contains ! Generic loop iterator integer :: i, j, k + integer :: ierr !< Generic flag used to identify and report database errors + ! Silo-HDF5 Database Format if (format == 1) then @@ -1106,7 +1112,8 @@ contains logical :: lg_bub_file, file_exist integer, dimension(2) :: gsizes, lsizes, start_idx_part - integer :: ifile, ierr, tot_data + integer :: ifile, tot_data + integer :: ierr !< Generic flag used to identify and report MPI errors integer :: i write (file_loc, '(A,I0,A)') 'lag_bubbles_mpi_io_', t_step, '.dat' @@ -1392,6 +1399,8 @@ contains ! domain, because one associated with the entire domain is ! not generated. + integer :: ierr !< Generic flag used to identify and report database errors + ! Silo-HDF5 database format if (format == 1) then ierr = DBCLOSE(dbfile) diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 809bdb60a3..d7bdd89ab5 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -197,8 +197,6 @@ module m_global_parameters integer :: mpi_info_int !> @} - integer, private :: ierr - type(physical_parameters), dimension(num_fluids_max) :: fluid_pp !< !! Database of the physical parameters of each of the fluids that is present !! in the flow. These include the stiffened gas equation of state parameters, @@ -904,6 +902,10 @@ contains !> Subroutine to initialize parallel infrastructure impure subroutine s_initialize_parallel_io +#ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors +#endif + num_dims = 1 + min(1, n) + min(1, p) if (mhd) then diff --git a/src/post_process/m_mpi_proxy.fpp b/src/post_process/m_mpi_proxy.fpp index 87f7f80d9f..4bd4987ca9 100644 --- a/src/post_process/m_mpi_proxy.fpp +++ b/src/post_process/m_mpi_proxy.fpp @@ -31,11 +31,6 @@ module m_mpi_proxy integer, allocatable, dimension(:) :: displs !> @} - !> @name Generic flags used to identify and report MPI errors - !> @{ - integer, private :: ierr - !> @} - contains !> Computation of parameters, allocation procedures, and/or @@ -45,6 +40,7 @@ contains #ifdef MFC_MPI integer :: i !< Generic loop iterator + integer :: ierr !< Generic flag used to identify and report MPI errors ! Allocating and configuring the receive counts and the displacement ! vector variables used in variable-gather communication procedures. @@ -85,6 +81,7 @@ contains #ifdef MFC_MPI integer :: i !< Generic loop iterator + integer :: ierr !< Generic flag used to identify and report MPI errors ! Logistics call MPI_BCAST(case_dir, len(case_dir), MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr) @@ -150,6 +147,7 @@ contains real(wp), dimension(1:, 0:), intent(INOUT) :: spatial_extents #ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors ! Simulation is 3D if (p > 0) then @@ -267,6 +265,7 @@ contains impure subroutine s_mpi_defragment_1d_grid_variable #ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors ! Silo-HDF5 database format if (format == 1) then @@ -306,6 +305,7 @@ contains real(wp), dimension(1:2, 0:num_procs - 1), intent(inout) :: data_extents #ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors ! Minimum flow variable extent call MPI_GATHERV(minval(q_sf), 1, mpi_p, & @@ -333,6 +333,7 @@ contains real(wp), dimension(0:m), intent(inout) :: q_root_sf #ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors ! Gathering the sub-domain flow variable data from all the processes ! and putting it back together for the entire computational domain diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index 9d3832f680..9a5eb51391 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -188,8 +188,6 @@ module m_global_parameters #endif - integer, private :: ierr - ! Initial Condition Parameters integer :: num_patches !< Number of patches composing initial condition @@ -931,6 +929,10 @@ contains impure subroutine s_initialize_parallel_io +#ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors +#endif + num_dims = 1 + min(1, n) + min(1, p) if (mhd) then diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp index 82295cd4ed..1a0e691aa4 100644 --- a/src/pre_process/m_mpi_proxy.fpp +++ b/src/pre_process/m_mpi_proxy.fpp @@ -23,9 +23,6 @@ module m_mpi_proxy implicit none - integer, private :: ierr !< - !! Generic flag used to identify and report MPI errors - contains !> Since only processor with rank 0 is in charge of reading !! and checking the consistency of the user provided inputs, @@ -38,6 +35,8 @@ contains ! Generic loop iterator integer :: i + ! Generic flag used to identify and report MPI errors + integer :: ierr ! Logistics call MPI_BCAST(case_dir, len(case_dir), MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr) diff --git a/src/pre_process/m_patches.fpp b/src/pre_process/m_patches.fpp index 0d512c5d9f..db09af9d1f 100644 --- a/src/pre_process/m_patches.fpp +++ b/src/pre_process/m_patches.fpp @@ -380,12 +380,12 @@ contains !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables !! @param ib True if this patch is an immersed boundary - subroutine s_circle(patch_id, patch_id_fp, q_prim_vf, ib) + subroutine s_circle(patch_id, patch_id_fp, q_prim_vf, ib_flag) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - logical, optional, intent(in) :: ib + logical, optional, intent(in) :: ib_flag real(wp) :: radius @@ -396,7 +396,7 @@ contains ! Transferring the circular patch's radius, centroid, smearing patch ! identity and smearing coefficient information - if (present(ib)) then + if (present(ib_flag)) then x_centroid = patch_ib(patch_id)%x_centroid y_centroid = patch_ib(patch_id)%y_centroid radius = patch_ib(patch_id)%radius @@ -421,7 +421,7 @@ contains do j = 0, n do i = 0, m - if (.not. present(ib) .and. patch_icpp(patch_id)%smoothen) then + if (.not. present(ib_flag) .and. patch_icpp(patch_id)%smoothen) then eta = tanh(smooth_coeff/min(dx, dy)* & (sqrt((x_cc(i) - x_centroid)**2 & @@ -430,8 +430,8 @@ contains end if - if (present(ib) .and. ((x_cc(i) - x_centroid)**2 & - + (y_cc(j) - y_centroid)**2 <= radius**2)) & + if (present(ib_flag) .and. ((x_cc(i) - x_centroid)**2 & + + (y_cc(j) - y_centroid)**2 <= radius**2)) & then patch_id_fp(i, j, 0) = patch_id @@ -441,7 +441,7 @@ contains .and. & patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) & .or. & - (.not. present(ib) .and. patch_id_fp(i, j, 0) == smooth_patch_id)) & + (.not. present(ib_flag) .and. patch_id_fp(i, j, 0) == smooth_patch_id)) & then call s_assign_patch_primitive_variables(patch_id, i, j, 0, & @@ -464,28 +464,29 @@ contains !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables !! @param ib True if this patch is an immersed boundary - subroutine s_airfoil(patch_id, patch_id_fp, q_prim_vf, ib) + subroutine s_airfoil(patch_id, patch_id_fp, q_prim_vf, ib_flag) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - logical, optional, intent(in) :: ib + logical, optional, intent(in) :: ib_flag - real(wp) :: x0, y0, f, x_act, y_act, ca, pa, ma, ta, theta, xa, yt, xu, yu, xl, yl, xc, yc, dycdxc, sin_c, cos_c + real(wp) :: x0, y0, f, x_act, y_act, ca_in, pa, ma, ta, theta + real(wp) :: xa, yt, xu, yu, xl, yl, xc, yc, dycdxc, sin_c, cos_c integer :: i, j, k integer :: Np1, Np2 - if (.not. present(ib)) return + if (.not. present(ib_flag)) return x0 = patch_ib(patch_id)%x_centroid y0 = patch_ib(patch_id)%y_centroid - ca = patch_ib(patch_id)%c + ca_in = patch_ib(patch_id)%c pa = patch_ib(patch_id)%p ma = patch_ib(patch_id)%m ta = patch_ib(patch_id)%t theta = pi*patch_ib(patch_id)%theta/180._wp - Np1 = int((pa*ca/dx)*20) - Np2 = int(((ca - pa*ca)/dx)*20) + Np1 = int((pa*ca_in/dx)*20) + Np2 = int(((ca_in - pa*ca_in)/dx)*20) Np = Np1 + Np2 + 1 allocate (airfoil_grid_u(1:Np)) @@ -501,13 +502,13 @@ contains do i = 1, Np1 + Np2 - 1 if (i <= Np1) then - xc = x0 + i*(pa*ca/Np1) - xa = (xc - x0)/ca + xc = x0 + i*(pa*ca_in/Np1) + xa = (xc - x0)/ca_in yc = (ma/pa**2)*(2*pa*xa - xa**2) dycdxc = (2*ma/pa**2)*(pa - xa) else - xc = x0 + pa*ca + (i - Np1)*((ca - pa*ca)/Np2) - xa = (xc - x0)/ca + xc = x0 + pa*ca_in + (i - Np1)*((ca_in - pa*ca_in)/Np2) + xa = (xc - x0)/ca_in yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2) dycdxc = (2*ma/(1 - pa)**2)*(pa - xa) end if @@ -522,11 +523,11 @@ contains xl = xa + yt*sin_c yl = yc - yt*cos_c - xu = xu*ca + x0 - yu = yu*ca + y0 + xu = xu*ca_in + x0 + yu = yu*ca_in + y0 - xl = xl*ca + x0 - yl = yl*ca + y0 + xl = xl*ca_in + x0 + yl = yl*ca_in + y0 airfoil_grid_u(i + 1)%x = xu airfoil_grid_u(i + 1)%y = yu @@ -536,10 +537,10 @@ contains end do - airfoil_grid_u(Np)%x = x0 + ca + airfoil_grid_u(Np)%x = x0 + ca_in airfoil_grid_u(Np)%y = y0 - airfoil_grid_l(Np)%x = x0 + ca + airfoil_grid_l(Np)%x = x0 + ca_in airfoil_grid_l(Np)%y = y0 do j = 0, n @@ -553,8 +554,8 @@ contains y_act = y_cc(j) end if - if (x_act >= x0 .and. x_act <= x0 + ca) then - xa = (x_act - x0)/ca + if (x_act >= x0 .and. x_act <= x0 + ca_in) then + xa = (x_act - x0)/ca_in if (xa <= pa) then yc = (ma/pa**2)*(2*pa*xa - xa**2) dycdxc = (2*ma/pa**2)*(pa - xa) @@ -626,30 +627,30 @@ contains !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables !! @param ib True if this patch is an immersed boundary - subroutine s_3D_airfoil(patch_id, patch_id_fp, q_prim_vf, ib) + subroutine s_3D_airfoil(patch_id, patch_id_fp, q_prim_vf, ib_flag) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - logical, optional, intent(in) :: ib + logical, optional, intent(in) :: ib_flag - real(wp) :: x0, y0, z0, lz, z_max, z_min, f, x_act, y_act, ca, pa, ma, ta, theta, xa, yt, xu, yu, xl, yl, xc, yc, dycdxc, sin_c, cos_c + real(wp) :: x0, y0, z0, lz, z_max, z_min, f, x_act, y_act, ca_in, pa, ma, ta, theta, xa, yt, xu, yu, xl, yl, xc, yc, dycdxc, sin_c, cos_c integer :: i, j, k, l integer :: Np1, Np2 - if (.not. present(ib)) return + if (.not. present(ib_flag)) return x0 = patch_ib(patch_id)%x_centroid y0 = patch_ib(patch_id)%y_centroid z0 = patch_ib(patch_id)%z_centroid lz = patch_ib(patch_id)%length_z - ca = patch_ib(patch_id)%c + ca_in = patch_ib(patch_id)%c pa = patch_ib(patch_id)%p ma = patch_ib(patch_id)%m ta = patch_ib(patch_id)%t theta = pi*patch_ib(patch_id)%theta/180._wp - Np1 = int((pa*ca/dx)*20) - Np2 = int(((ca - pa*ca)/dx)*20) + Np1 = int((pa*ca_in/dx)*20) + Np2 = int(((ca_in - pa*ca_in)/dx)*20) Np = Np1 + Np2 + 1 allocate (airfoil_grid_u(1:Np)) @@ -668,13 +669,13 @@ contains do i = 1, Np1 + Np2 - 1 if (i <= Np1) then - xc = x0 + i*(pa*ca/Np1) - xa = (xc - x0)/ca + xc = x0 + i*(pa*ca_in/Np1) + xa = (xc - x0)/ca_in yc = (ma/pa**2)*(2*pa*xa - xa**2) dycdxc = (2*ma/pa**2)*(pa - xa) else - xc = x0 + pa*ca + (i - Np1)*((ca - pa*ca)/Np2) - xa = (xc - x0)/ca + xc = x0 + pa*ca_in + (i - Np1)*((ca_in - pa*ca_in)/Np2) + xa = (xc - x0)/ca_in yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2) dycdxc = (2*ma/(1 - pa)**2)*(pa - xa) end if @@ -689,11 +690,11 @@ contains xl = xa + yt*sin_c yl = yc - yt*cos_c - xu = xu*ca + x0 - yu = yu*ca + y0 + xu = xu*ca_in + x0 + yu = yu*ca_in + y0 - xl = xl*ca + x0 - yl = yl*ca + y0 + xl = xl*ca_in + x0 + yl = yl*ca_in + y0 airfoil_grid_u(i + 1)%x = xu airfoil_grid_u(i + 1)%y = yu @@ -703,10 +704,10 @@ contains end do - airfoil_grid_u(Np)%x = x0 + ca + airfoil_grid_u(Np)%x = x0 + ca_in airfoil_grid_u(Np)%y = y0 - airfoil_grid_l(Np)%x = x0 + ca + airfoil_grid_l(Np)%x = x0 + ca_in airfoil_grid_l(Np)%y = y0 do l = 0, p @@ -722,8 +723,8 @@ contains y_act = y_cc(j) end if - if (x_act >= x0 .and. x_act <= x0 + ca) then - xa = (x_act - x0)/ca + if (x_act >= x0 .and. x_act <= x0 + ca_in) then + xa = (x_act - x0)/ca_in if (xa <= pa) then yc = (ma/pa**2)*(2*pa*xa - xa**2) dycdxc = (2*ma/pa**2)*(pa - xa) @@ -1100,12 +1101,12 @@ contains !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables !! @param ib True if this patch is an immersed boundary - subroutine s_rectangle(patch_id, patch_id_fp, q_prim_vf, ib) + subroutine s_rectangle(patch_id, patch_id_fp, q_prim_vf, ib_flag) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - logical, optional, intent(in) :: ib !< True if this patch is an immersed boundary + logical, optional, intent(in) :: ib_flag !< True if this patch is an immersed boundary integer :: i, j, k !< generic loop iterators real(wp) :: pi_inf, gamma, lit_gamma !< Equation of state parameters @@ -1117,7 +1118,7 @@ contains lit_gamma = (1._wp + gamma)/gamma ! Transferring the rectangle's centroid and length information - if (present(ib)) then + if (present(ib_flag)) then x_centroid = patch_ib(patch_id)%x_centroid y_centroid = patch_ib(patch_id)%y_centroid length_x = patch_ib(patch_id)%length_x @@ -1152,7 +1153,7 @@ contains x_boundary%end >= x_cc(i) .and. & y_boundary%beg <= y_cc(j) .and. & y_boundary%end >= y_cc(j)) then - if (present(ib)) then + if (present(ib_flag)) then ! Updating the patch identities bookkeeping variable patch_id_fp(i, j, 0) = patch_id else @@ -1417,17 +1418,17 @@ contains real(wp) :: r, x_p, eps, phi real(wp), dimension(2:9) :: as, Ps - real(wp) :: radius, x_centroid, y_centroid, z_centroid, eta, smooth_coeff - logical :: non_axis_sym + real(wp) :: radius, x_centroid_local, y_centroid_local, z_centroid_local, eta_local, smooth_coeff_local + logical :: non_axis_sym_in integer :: i, j, k !< generic loop iterators ! Transferring the patch's centroid and radius information - x_centroid = patch_icpp(patch_id)%x_centroid - y_centroid = patch_icpp(patch_id)%y_centroid - z_centroid = patch_icpp(patch_id)%z_centroid + x_centroid_local = patch_icpp(patch_id)%x_centroid + y_centroid_local = patch_icpp(patch_id)%y_centroid + z_centroid_local = patch_icpp(patch_id)%z_centroid smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id - smooth_coeff = patch_icpp(patch_id)%smooth_coeff + smooth_coeff_local = patch_icpp(patch_id)%smooth_coeff radius = patch_icpp(patch_id)%radius as(2) = patch_icpp(patch_id)%a(2) as(3) = patch_icpp(patch_id)%a(3) @@ -1437,20 +1438,20 @@ contains as(7) = patch_icpp(patch_id)%a(7) as(8) = patch_icpp(patch_id)%a(8) as(9) = patch_icpp(patch_id)%a(9) - non_axis_sym = patch_icpp(patch_id)%non_axis_sym + non_axis_sym_in = patch_icpp(patch_id)%non_axis_sym ! Since the analytical patch does not allow for its boundaries to get ! smoothed out, the pseudo volume fraction is set to 1 to make sure ! that only the current patch contributes to the fluid state in the ! cells that this patch covers. - eta = 1._wp + eta_local = 1._wp eps = 1.e-32_wp ! Checking whether the patch covers a particular cell in the domain ! and verifying whether the current patch has permission to write to ! to that cell. If both queries check out, the primitive variables ! of the current patch are assigned to this cell. - if (p > 0 .and. .not. non_axis_sym) then + if (p > 0 .and. .not. non_axis_sym_in) then do k = 0, p do j = 0, n do i = 0, m @@ -1461,11 +1462,11 @@ contains cart_z = z_cc(k) end if - r = sqrt((x_cc(i) - x_centroid)**2 + (cart_y - y_centroid)**2 + (cart_z - z_centroid)**2) + eps - if (x_cc(i) - x_centroid <= 0) then - x_p = -1._wp*abs(x_cc(i) - x_centroid + eps)/r + r = sqrt((x_cc(i) - x_centroid_local)**2 + (cart_y - y_centroid_local)**2 + (cart_z - z_centroid_local)**2) + eps + if (x_cc(i) - x_centroid_local <= 0) then + x_p = -1._wp*abs(x_cc(i) - x_centroid_local + eps)/r else - x_p = abs(x_cc(i) - x_centroid + eps)/r + x_p = abs(x_cc(i) - x_centroid_local + eps)/r end if Ps(2) = unassociated_legendre(x_p, 2) @@ -1474,7 +1475,7 @@ contains Ps(5) = unassociated_legendre(x_p, 5) Ps(6) = unassociated_legendre(x_p, 6) Ps(7) = unassociated_legendre(x_p, 7) - if ((x_cc(i) - x_centroid >= 0 & + if ((x_cc(i) - x_centroid_local >= 0 & .and. & r - as(2)*Ps(2) - as(3)*Ps(3) - as(4)*Ps(4) - as(5)*Ps(5) - as(6)*Ps(6) - as(7)*Ps(7) <= radius & .and. & @@ -1482,13 +1483,13 @@ contains (patch_id_fp(i, j, k) == smooth_patch_id)) & then if (patch_icpp(patch_id)%smoothen) then - eta = tanh(smooth_coeff/min(dx, dy, dz)* & - ((r - as(2)*Ps(2) - as(3)*Ps(3) - as(4)*Ps(4) - as(5)*Ps(5) - as(6)*Ps(6) - as(7)*Ps(7)) & - - radius))*(-0.5_wp) + 0.5_wp + eta_local = tanh(smooth_coeff_local/min(dx, dy, dz)* & + ((r - as(2)*Ps(2) - as(3)*Ps(3) - as(4)*Ps(4) - as(5)*Ps(5) - as(6)*Ps(6) - as(7)*Ps(7)) & + - radius))*(-0.5_wp) + 0.5_wp end if call s_assign_patch_primitive_variables(patch_id, i, j, k, & - eta, q_prim_vf, patch_id_fp) + eta_local, q_prim_vf, patch_id_fp) end if end do @@ -1499,9 +1500,9 @@ contains do j = 0, n do i = 0, m - if (non_axis_sym) then - phi = atan(((y_cc(j) - y_centroid) + eps)/((x_cc(i) - x_centroid) + eps)) - r = sqrt((x_cc(i) - x_centroid)**2._wp + (y_cc(j) - y_centroid)**2._wp) + eps + if (non_axis_sym_in) then + phi = atan(((y_cc(j) - y_centroid_local) + eps)/((x_cc(i) - x_centroid_local) + eps)) + r = sqrt((x_cc(i) - x_centroid_local)**2._wp + (y_cc(j) - y_centroid_local)**2._wp) + eps x_p = (eps)/r Ps(2) = spherical_harmonic_func(x_p, phi, 2, 2) Ps(3) = spherical_harmonic_func(x_p, phi, 3, 3) @@ -1512,8 +1513,8 @@ contains Ps(8) = spherical_harmonic_func(x_p, phi, 8, 8) Ps(9) = spherical_harmonic_func(x_p, phi, 9, 9) else - r = sqrt((x_cc(i) - x_centroid)**2._wp + (y_cc(j) - y_centroid)**2._wp) + eps - x_p = abs(x_cc(i) - x_centroid + eps)/r + r = sqrt((x_cc(i) - x_centroid_local)**2._wp + (y_cc(j) - y_centroid_local)**2._wp) + eps + x_p = abs(x_cc(i) - x_centroid_local + eps)/r Ps(2) = unassociated_legendre(x_p, 2) Ps(3) = unassociated_legendre(x_p, 3) Ps(4) = unassociated_legendre(x_p, 4) @@ -1524,22 +1525,22 @@ contains Ps(9) = unassociated_legendre(x_p, 9) end if - if (x_cc(i) - x_centroid >= 0 & + if (x_cc(i) - x_centroid_local >= 0 & .and. & r - as(2)*Ps(2) - as(3)*Ps(3) - as(4)*Ps(4) - as(5)*Ps(5) - as(6)*Ps(6) - as(7)*Ps(7) - as(8)*Ps(8) - as(9)*Ps(9) <= radius .and. & patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) & then call s_assign_patch_primitive_variables(patch_id, i, j, 0, & - eta, q_prim_vf, patch_id_fp) + eta_local, q_prim_vf, patch_id_fp) - elseif (x_cc(i) - x_centroid < 0 & + elseif (x_cc(i) - x_centroid_local < 0 & .and. & r - as(2)*Ps(2) + as(3)*Ps(3) - as(4)*Ps(4) + as(5)*Ps(5) - as(6)*Ps(6) + as(7)*Ps(7) - as(8)*Ps(8) + as(9)*Ps(9) <= radius & .and. & patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) & then call s_assign_patch_primitive_variables(patch_id, i, j, 0, & - eta, q_prim_vf, patch_id_fp) + eta_local, q_prim_vf, patch_id_fp) end if end do @@ -1557,12 +1558,12 @@ contains !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables !! @param ib True if this patch is an immersed boundary - subroutine s_sphere(patch_id, patch_id_fp, q_prim_vf, ib) + subroutine s_sphere(patch_id, patch_id_fp, q_prim_vf, ib_flag) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - logical, optional, intent(in) :: ib !< True if this patch is an immersed boundary + logical, optional, intent(in) :: ib_flag !< True if this patch is an immersed boundary ! Generic loop iterators integer :: i, j, k @@ -1575,7 +1576,7 @@ contains ! Transferring spherical patch's radius, centroid, smoothing patch ! identity and smoothing coefficient information - if (present(ib)) then + if (present(ib_flag)) then x_centroid = patch_ib(patch_id)%x_centroid y_centroid = patch_ib(patch_id)%y_centroid z_centroid = patch_ib(patch_id)%z_centroid @@ -1609,7 +1610,7 @@ contains cart_z = z_cc(k) end if - if (.not. present(ib) .and. patch_icpp(patch_id)%smoothen) then + if (.not. present(ib_flag) .and. patch_icpp(patch_id)%smoothen) then eta = tanh(smooth_coeff/min(dx, dy, dz)* & (sqrt((x_cc(i) - x_centroid)**2 & + (cart_y - y_centroid)**2 & @@ -1617,7 +1618,7 @@ contains - radius))*(-0.5_wp) + 0.5_wp end if - if (present(ib)) then + if (present(ib_flag)) then ! Updating the patch identities bookkeeping variable if (((x_cc(i) - x_centroid)**2 & + (cart_y - y_centroid)**2 & @@ -1659,10 +1660,10 @@ contains !! @param patch_id is the patch identifier !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables - subroutine s_cuboid(patch_id, patch_id_fp, q_prim_vf, ib) + subroutine s_cuboid(patch_id, patch_id_fp, q_prim_vf, ib_flag) integer, intent(in) :: patch_id - logical, optional, intent(in) :: ib + logical, optional, intent(in) :: ib_flag integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf @@ -1671,7 +1672,7 @@ contains @:Hardcoded3DVariables() ! Transferring the cuboid's centroid and length information - if (present(ib)) then + if (present(ib_flag)) then x_centroid = patch_ib(patch_id)%x_centroid y_centroid = patch_ib(patch_id)%y_centroid z_centroid = patch_ib(patch_id)%z_centroid @@ -1724,7 +1725,7 @@ contains z_boundary%beg <= cart_z .and. & z_boundary%end >= cart_z) then - if (present(ib)) then + if (present(ib_flag)) then ! Updating the patch identities bookkeeping variable patch_id_fp(i, j, k) = patch_id else @@ -1763,12 +1764,12 @@ contains !! @param patch_id_fp Array to track patch ids !! @param q_prim_vf Array of primitive variables !! @param ib True if this patch is an immersed boundary - subroutine s_cylinder(patch_id, patch_id_fp, q_prim_vf, ib) + subroutine s_cylinder(patch_id, patch_id_fp, q_prim_vf, ib_flag) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp type(scalar_field), dimension(1:sys_size), intent(inout) :: q_prim_vf - logical, optional, intent(in) :: ib !< True if this patch is an immersed boundary + logical, optional, intent(in) :: ib_flag !< True if this patch is an immersed boundary integer :: i, j, k !< Generic loop iterators real(wp) :: radius @@ -1778,7 +1779,7 @@ contains ! Transferring the cylindrical patch's centroid, length, radius, ! smoothing patch identity and smoothing coefficient information - if (present(ib)) then + if (present(ib_flag)) then x_centroid = patch_ib(patch_id)%x_centroid y_centroid = patch_ib(patch_id)%y_centroid z_centroid = patch_ib(patch_id)%z_centroid @@ -1827,7 +1828,7 @@ contains cart_z = z_cc(k) end if - if (.not. present(ib) .and. patch_icpp(patch_id)%smoothen) then + if (.not. present(ib_flag) .and. patch_icpp(patch_id)%smoothen) then if (.not. f_is_default(length_x)) then eta = tanh(smooth_coeff/min(dy, dz)* & (sqrt((cart_y - y_centroid)**2 & @@ -1846,7 +1847,7 @@ contains end if end if - if (present(ib)) then + if (present(ib_flag)) then if (((.not. f_is_default(length_x) .and. & (cart_y - y_centroid)**2 & + (cart_z - z_centroid)**2 <= radius**2 .and. & @@ -2004,7 +2005,7 @@ contains !! @param ib True if this patch is an immersed boundary !! @param STL_levelset STL levelset !! @param STL_levelset_norm STL levelset normals - subroutine s_model(patch_id, patch_id_fp, q_prim_vf, ib, STL_levelset, STL_levelset_norm) + subroutine s_model(patch_id, patch_id_fp, q_prim_vf, ib_flag, STL_levelset, STL_levelset_norm) integer, intent(in) :: patch_id integer, dimension(0:m, 0:n, 0:p), intent(inout) :: patch_id_fp @@ -2013,7 +2014,7 @@ contains ! Variables for IBM+STL type(levelset_field), optional, intent(inout) :: STL_levelset !< Levelset determined by models type(levelset_norm_field), optional, intent(inout) :: STL_levelset_norm !< Levelset_norm determined by models - logical, optional, intent(in) :: ib !< True if this patch is an immersed boundary + logical, optional, intent(in) :: ib_flag !< True if this patch is an immersed boundary real(wp) :: normals(1:3) !< Boundary normal buffer integer :: boundary_vertex_count, boundary_edge_count, total_vertices !< Boundary vertex real(wp), allocatable, dimension(:, :, :) :: boundary_v !< Boundary vertex buffer @@ -2036,13 +2037,13 @@ contains real(wp), dimension(1:4, 1:4) :: transform, transform_n - if (present(ib) .and. proc_rank == 0) then + if (present(ib_flag) .and. proc_rank == 0) then print *, " * Reading model: "//trim(patch_ib(patch_id)%model_filepath) else if (proc_rank == 0) then print *, " * Reading model: "//trim(patch_icpp(patch_id)%model_filepath) end if - if (present(ib)) then + if (present(ib_flag)) then model = f_model_read(patch_ib(patch_id)%model_filepath) params%scale(:) = patch_ib(patch_id)%model_scale(:) params%translate(:) = patch_ib(patch_id)%model_translate(:) @@ -2152,13 +2153,13 @@ contains point = f_convert_cyl_to_cart(point) end if - if (present(ib)) then + if (present(ib_flag)) then eta = f_model_is_inside(model, point, (/dx, dy, dz/), patch_ib(patch_id)%model_spc) else eta = f_model_is_inside(model, point, (/dx, dy, dz/), patch_icpp(patch_id)%model_spc) end if - if (present(ib)) then + if (present(ib_flag)) then ! Reading STL boundary vertices and compute the levelset and levelset_norm if (eta > patch_ib(patch_id)%model_threshold) then patch_id_fp(i, j, k) = patch_id diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index 3a02f2f69e..1779b464cf 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -82,16 +82,16 @@ module m_start_up !! @param q_cons_vf Conservative variables !! @param ib_markers track if a cell is within the immersed boundary - impure subroutine s_read_abstract_ic_data_files(q_cons_vf, ib_markers) + impure subroutine s_read_abstract_ic_data_files(q_cons_vf_in, ib_markers_in) import :: scalar_field, integer_field, sys_size, pres_field type(scalar_field), & dimension(sys_size), & - intent(inout) :: q_cons_vf + intent(inout) :: q_cons_vf_in type(integer_field), & - intent(inout) :: ib_markers + intent(inout) :: ib_markers_in end subroutine s_read_abstract_ic_data_files @@ -409,14 +409,14 @@ contains !! all new initial condition. !! @param q_cons_vf Conservative variables !! @param ib_markers track if a cell is within the immersed boundary - impure subroutine s_read_serial_ic_data_files(q_cons_vf, ib_markers) + impure subroutine s_read_serial_ic_data_files(q_cons_vf_in, ib_markers_in) type(scalar_field), & dimension(sys_size), & - intent(inout) :: q_cons_vf + intent(inout) :: q_cons_vf_in type(integer_field), & - intent(inout) :: ib_markers + intent(inout) :: ib_markers_in character(LEN=len_trim(case_dir) + 3*name_len) :: file_loc !< ! Generic string used to store the address of a particular file @@ -446,7 +446,7 @@ contains if (file_check) then open (1, FILE=trim(file_loc), FORM='unformatted', & STATUS='old', ACTION='read') - read (1) q_cons_vf(i)%sf + read (1) q_cons_vf_in(i)%sf close (1) else call s_mpi_abort('File q_cons_vf'//trim(file_num)// & @@ -517,7 +517,7 @@ contains if (file_check) then open (1, FILE=trim(file_loc), FORM='unformatted', & STATUS='old', ACTION='read') - read (1) ib_markers%sf(0:m, 0:n, 0:p) + read (1) ib_markers_in%sf(0:m, 0:n, 0:p) close (1) else call s_mpi_abort('File ib.dat is missing in ' & @@ -645,14 +645,14 @@ contains !! all new initial condition. !! @param q_cons_vf Conservative variables !! @param ib_markers track if a cell is within the immersed boundary - impure subroutine s_read_parallel_ic_data_files(q_cons_vf, ib_markers) + impure subroutine s_read_parallel_ic_data_files(q_cons_vf_in, ib_markers_in) type(scalar_field), & dimension(sys_size), & - intent(inout) :: q_cons_vf + intent(inout) :: q_cons_vf_in type(integer_field), & - intent(inout) :: ib_markers + intent(inout) :: ib_markers_in #ifdef MFC_MPI @@ -683,9 +683,9 @@ contains ! Initialize MPI data I/O if (ib) then - call s_initialize_mpi_data(q_cons_vf, ib_markers, levelset, levelset_norm) + call s_initialize_mpi_data(q_cons_vf_in, ib_markers_in, levelset, levelset_norm) else - call s_initialize_mpi_data(q_cons_vf) + call s_initialize_mpi_data(q_cons_vf_in) end if ! Size of local arrays diff --git a/src/simulation/include/inline_capillary.fpp b/src/simulation/include/inline_capillary.fpp index c36c617a27..417ff8977d 100644 --- a/src/simulation/include/inline_capillary.fpp +++ b/src/simulation/include/inline_capillary.fpp @@ -1,4 +1,4 @@ -#:def compute_capilary_stress_tensor() +#:def compute_capillary_stress_tensor() Omega(1, 1) = -sigma*(w2*w2 + w3*w3)/normW @@ -19,4 +19,4 @@ end if -#:enddef compute_capilary_stress_tensor +#:enddef compute_capillary_stress_tensor diff --git a/src/simulation/m_bubbles.fpp b/src/simulation/m_bubbles.fpp index 6ee19c210c..0ec758dc22 100644 --- a/src/simulation/m_bubbles.fpp +++ b/src/simulation/m_bubbles.fpp @@ -316,23 +316,23 @@ contains !> Subroutine that computes bubble wall properties for vapor bubbles !! @param pb Internal bubble pressure !! @param iR0 Current bubble size index - pure elemental subroutine s_bwproperty(pb, iR0, chi_vw, k_mw, rho_mw) + pure elemental subroutine s_bwproperty(pb_in, iR0, chi_vw_out, k_mw_out, rho_mw_out) $:GPU_ROUTINE(parallelism='[seq]') - real(wp), intent(in) :: pb + real(wp), intent(in) :: pb_in integer, intent(in) :: iR0 - real(wp), intent(out) :: chi_vw !< Bubble wall properties (Ando 2010) - real(wp), intent(out) :: k_mw !< Bubble wall properties (Ando 2010) - real(wp), intent(out) :: rho_mw !< Bubble wall properties (Ando 2010) + real(wp), intent(out) :: chi_vw_out !< Bubble wall properties (Ando 2010) + real(wp), intent(out) :: k_mw_out !< Bubble wall properties (Ando 2010) + real(wp), intent(out) :: rho_mw_out !< Bubble wall properties (Ando 2010) real(wp) :: x_vw ! mass fraction of vapor - chi_vw = 1._wp/(1._wp + R_v/R_n*(pb/pv - 1._wp)) + chi_vw_out = 1._wp/(1._wp + R_v/R_n*(pb_in/pv - 1._wp)) ! mole fraction of vapor & thermal conductivity of gas mixture - x_vw = M_n*chi_vw/(M_v + (M_n - M_v)*chi_vw) - k_mw = x_vw*k_v(iR0)/(x_vw + (1._wp - x_vw)*phi_vn) & - + (1._wp - x_vw)*k_n(iR0)/(x_vw*phi_nv + 1._wp - x_vw) + x_vw = M_n*chi_vw_out/(M_v + (M_n - M_v)*chi_vw_out) + k_mw_out = x_vw*k_v(iR0)/(x_vw + (1._wp - x_vw)*phi_vn) & + + (1._wp - x_vw)*k_n(iR0)/(x_vw*phi_nv + 1._wp - x_vw) ! gas mixture density - rho_mw = pv/(chi_vw*R_v*Tw) + rho_mw_out = pv/(chi_vw_out*R_v*Tw) end subroutine s_bwproperty diff --git a/src/simulation/m_bubbles_EE.fpp b/src/simulation/m_bubbles_EE.fpp index 5fe6d38136..b43a89e2e5 100644 --- a/src/simulation/m_bubbles_EE.fpp +++ b/src/simulation/m_bubbles_EE.fpp @@ -92,11 +92,11 @@ contains end subroutine s_comp_alpha_from_n - pure subroutine s_compute_bubbles_EE_rhs(idir, q_prim_vf, divu) + pure subroutine s_compute_bubbles_EE_rhs(idir, q_prim_vf, divu_in) integer, intent(in) :: idir type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf - type(scalar_field), intent(inout) :: divu !< matrix for div(u) + type(scalar_field), intent(inout) :: divu_in !< matrix for div(u) integer :: j, k, l @@ -107,8 +107,8 @@ contains do l = 0, p do k = 0, n do j = 0, m - divu%sf(j, k, l) = 0._wp - divu%sf(j, k, l) = & + divu_in%sf(j, k, l) = 0._wp + divu_in%sf(j, k, l) = & 5.e-1_wp/dx(j)*(q_prim_vf(contxe + idir)%sf(j + 1, k, l) - & q_prim_vf(contxe + idir)%sf(j - 1, k, l)) @@ -123,9 +123,9 @@ contains do l = 0, p do k = 0, n do j = 0, m - divu%sf(j, k, l) = divu%sf(j, k, l) + & - 5.e-1_wp/dy(k)*(q_prim_vf(contxe + idir)%sf(j, k + 1, l) - & - q_prim_vf(contxe + idir)%sf(j, k - 1, l)) + divu_in%sf(j, k, l) = divu_in%sf(j, k, l) + & + 5.e-1_wp/dy(k)*(q_prim_vf(contxe + idir)%sf(j, k + 1, l) - & + q_prim_vf(contxe + idir)%sf(j, k - 1, l)) end do end do @@ -137,9 +137,9 @@ contains do l = 0, p do k = 0, n do j = 0, m - divu%sf(j, k, l) = divu%sf(j, k, l) + & - 5.e-1_wp/dz(l)*(q_prim_vf(contxe + idir)%sf(j, k, l + 1) - & - q_prim_vf(contxe + idir)%sf(j, k, l - 1)) + divu_in%sf(j, k, l) = divu_in%sf(j, k, l) + & + 5.e-1_wp/dz(l)*(q_prim_vf(contxe + idir)%sf(j, k, l + 1) - & + q_prim_vf(contxe + idir)%sf(j, k, l - 1)) end do end do @@ -153,13 +153,14 @@ contains !! that are needed for the bubble modeling !! @param q_prim_vf Primitive variables !! @param q_cons_vf Conservative variables - impure subroutine s_compute_bubble_EE_source(q_cons_vf, q_prim_vf, rhs_vf) + impure subroutine s_compute_bubble_EE_source(q_cons_vf, q_prim_vf, rhs_vf, divu_in) type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf + type(scalar_field), intent(in) :: divu_in !< matrix for div(u) real(wp) :: rddot - real(wp) :: pb, mv, vflux, pbdot + real(wp) :: pb_local, mv_local, vflux, pbdot real(wp) :: n_tait, B_tait real(wp), dimension(nb) :: Rtmp, Vtmp real(wp) :: myR, myV, alf, myP, myRho, R2Vav, R3 @@ -270,16 +271,16 @@ contains myV = q_prim_vf(vs(q))%sf(j, k, l) if (.not. polytropic) then - pb = q_prim_vf(ps(q))%sf(j, k, l) - mv = q_prim_vf(ms(q))%sf(j, k, l) - call s_bwproperty(pb, q, chi_vw, k_mw, rho_mw) - call s_vflux(myR, myV, pb, mv, q, vflux) - pbdot = f_bpres_dot(vflux, myR, myV, pb, mv, q) + pb_local = q_prim_vf(ps(q))%sf(j, k, l) + mv_local = q_prim_vf(ms(q))%sf(j, k, l) + call s_bwproperty(pb_local, q, chi_vw, k_mw, rho_mw) + call s_vflux(myR, myV, pb_local, mv_local, q, vflux) + pbdot = f_bpres_dot(vflux, myR, myV, pb_local, mv_local, q) bub_p_src(j, k, l, q) = nbub*pbdot bub_m_src(j, k, l, q) = nbub*vflux*4._wp*pi*(myR**2._wp) else - pb = 0._wp; mv = 0._wp; vflux = 0._wp; pbdot = 0._wp + pb_local = 0._wp; mv_local = 0._wp; vflux = 0._wp; pbdot = 0._wp end if ! Adaptive time stepping @@ -288,8 +289,8 @@ contains if (adap_dt) then call s_advance_step(myRho, myP, myR, myV, R0(q), & - pb, pbdot, alf, n_tait, B_tait, & - bub_adv_src(j, k, l), divu%sf(j, k, l), & + pb_local, pbdot, alf, n_tait, B_tait, & + bub_adv_src(j, k, l), divu_in%sf(j, k, l), & dmBub_id, dmMass_v, dmMass_n, dmBeta_c, & dmBeta_t, dmCson, adap_dt_stop) @@ -298,8 +299,8 @@ contains else rddot = f_rddot(myRho, myP, myR, myV, R0(q), & - pb, pbdot, alf, n_tait, B_tait, & - bub_adv_src(j, k, l), divu%sf(j, k, l), & + pb_local, pbdot, alf, n_tait, B_tait, & + bub_adv_src(j, k, l), divu_in%sf(j, k, l), & dmCson) bub_v_src(j, k, l, q) = nbub*rddot bub_r_src(j, k, l, q) = q_cons_vf(vs(q))%sf(j, k, l) diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp index 20e7515c83..bacd19497d 100644 --- a/src/simulation/m_bubbles_EL.fpp +++ b/src/simulation/m_bubbles_EL.fpp @@ -293,7 +293,7 @@ contains integer :: i real(wp) :: pliq, volparticle, concvap, totalmass, kparticle, cpparticle - real(wp) :: omegaN, PeG, PeT, rhol, pcrit, qv, gamma, pi_inf, dynP + real(wp) :: omegaN_local, PeG, PeT, rhol, pcrit, qv, gamma, pi_inf, dynP integer, dimension(3) :: cell real(wp), dimension(2) :: Re real(wp) :: massflag, heatflag, Re_trans, Im_trans @@ -373,21 +373,21 @@ contains ! Bubble natural frequency concvap = gas_mv(bub_id, 1)/(gas_mv(bub_id, 1) + gas_mg(bub_id)) - omegaN = (3._wp*(gas_p(bub_id, 1) - pv*(massflag)) + 4._wp*(1._wp/Web)/bub_R0(bub_id))/rhol + omegaN_local = (3._wp*(gas_p(bub_id, 1) - pv*(massflag)) + 4._wp*(1._wp/Web)/bub_R0(bub_id))/rhol if (pv*(massflag) > gas_p(bub_id, 1)) then call s_mpi_abort("Lagrange bubble initially located in a region with pressure below the vapor pressure.") end if - omegaN = sqrt(omegaN/bub_R0(bub_id)**2._wp) + omegaN_local = sqrt(omegaN_local/bub_R0(bub_id)**2._wp) cpparticle = concvap*cp_v + (1._wp - concvap)*cp_n kparticle = concvap*k_vl + (1._wp - concvap)*k_nl ! Mass and heat transfer coefficients (based on Preston 2007) - PeT = totalmass/volparticle*cpparticle*bub_R0(bub_id)**2._wp*omegaN/kparticle + PeT = totalmass/volparticle*cpparticle*bub_R0(bub_id)**2._wp*omegaN_local/kparticle call s_transcoeff(1._wp, PeT, Re_trans, Im_trans) gas_betaT(bub_id) = Re_trans*(heatflag)*kparticle - PeG = bub_R0(bub_id)**2._wp*omegaN/lag_params%diffcoefvap + PeG = bub_R0(bub_id)**2._wp*omegaN_local/lag_params%diffcoefvap call s_transcoeff(1._wp, PeG, Re_trans, Im_trans) gas_betaC(bub_id) = Re_trans*(massflag)*lag_params%diffcoefvap diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index d655f1bbf6..da6cba12ac 100644 --- a/src/simulation/m_cbc.fpp +++ b/src/simulation/m_cbc.fpp @@ -649,7 +649,7 @@ contains real(wp), dimension(contxe) :: alpha_rho, dalpha_rho_ds, mf real(wp), dimension(2) :: Re_cbc real(wp), dimension(num_vels) :: vel, dvel_ds - real(wp), dimension(num_fluids) :: adv, dadv_ds + real(wp), dimension(num_fluids) :: adv_local, dadv_ds real(wp), dimension(sys_size) :: L real(wp), dimension(3) :: lambda @@ -769,7 +769,7 @@ contains end if ! FD2 or FD4 of RHS at j = 0 - $:GPU_PARALLEL_LOOP(collapse=2, private='[alpha_rho, vel, adv, & + $:GPU_PARALLEL_LOOP(collapse=2, private='[alpha_rho, vel, adv_local, & & mf, dvel_ds, dadv_ds, Re_cbc, dalpha_rho_ds,dvel_dt, & & dadv_dt, dalpha_rho_dt, L, lambda, Ys, dYs_dt, & & dYs_ds, h_k, Cp_i, Gamma_i, Xs]') @@ -797,13 +797,13 @@ contains $:GPU_LOOP(parallelism='[seq]') do i = 1, advxe - E_idx - adv(i) = q_prim_rs${XYZ}$_vf(0, k, r, E_idx + i) + adv_local(i) = q_prim_rs${XYZ}$_vf(0, k, r, E_idx + i) end do if (bubbles_euler) then - call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv, adv, alpha_rho, Re_cbc) + call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv, adv_local, alpha_rho, Re_cbc) else - call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, adv, alpha_rho, Re_cbc) + call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, adv_local, alpha_rho, Re_cbc) end if $:GPU_LOOP(parallelism='[seq]') @@ -841,7 +841,7 @@ contains H = (E + pres)/rho ! Compute mixture sound speed - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv, vel_K_sum, 0._wp, c) + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv_local, vel_K_sum, 0._wp, c) ! First-Order Spatial Derivatives of Primitive Variables @@ -934,7 +934,7 @@ contains end if $:GPU_LOOP(parallelism='[seq]') do i = E_idx, advxe - 1 - L(i) = c*Ma*(adv(i + 1 - E_idx) - alpha_in(i + 1 - E_idx, ${CBC_DIR}$))/Del_in(${CBC_DIR}$) + L(i) = c*Ma*(adv_local(i + 1 - E_idx) - alpha_in(i + 1 - E_idx, ${CBC_DIR}$))/Del_in(${CBC_DIR}$) end do L(advxe) = rho*c**2._wp*(1._wp + Ma)*(vel(dir_idx(1)) + vel_in(${CBC_DIR}$, dir_idx(1))*sign(1, cbc_loc))/Del_in(${CBC_DIR}$) + c*(1._wp + Ma)*(pres - pres_in(${CBC_DIR}$))/Del_in(${CBC_DIR}$) end if @@ -1003,7 +1003,7 @@ contains if (cyl_coord .and. cbc_dir == 2 .and. cbc_loc == 1) then $:GPU_LOOP(parallelism='[seq]') do i = 1, advxe - E_idx - dadv_dt(i) = -L(momxe + i) !+ adv(i) * vel(dir_idx(1))/y_cc(n) + dadv_dt(i) = -L(momxe + i) !+ adv_local(i) * vel(dir_idx(1))/y_cc(n) end do else $:GPU_LOOP(parallelism='[seq]') diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 63a7bad2c7..117429487b 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -1,4 +1,3 @@ - !! @file m_data_output.f90 !! @brief Contains module m_data_output @@ -985,10 +984,10 @@ contains !! @param t_step Current time-step !! @param q_com Center of mass information !! @param moments Higher moment information - impure subroutine s_write_com_files(t_step, c_mass) + impure subroutine s_write_com_files(t_step, c_mass_in) integer, intent(in) :: t_step - real(wp), dimension(num_fluids, 5), intent(in) :: c_mass + real(wp), dimension(num_fluids, 5), intent(in) :: c_mass_in integer :: i !< Generic loop iterator real(wp) :: nondim_time !< Non-dimensional time @@ -1004,28 +1003,28 @@ contains do i = 1, num_fluids ! Loop through fluids write (i + 120, '(6X,4F24.12)') & nondim_time, & - c_mass(i, 1), & - c_mass(i, 2), & - c_mass(i, 5) + c_mass_in(i, 1), & + c_mass_in(i, 2), & + c_mass_in(i, 5) end do elseif (p == 0) then ! 2D simulation do i = 1, num_fluids ! Loop through fluids write (i + 120, '(6X,5F24.12)') & nondim_time, & - c_mass(i, 1), & - c_mass(i, 2), & - c_mass(i, 3), & - c_mass(i, 5) + c_mass_in(i, 1), & + c_mass_in(i, 2), & + c_mass_in(i, 3), & + c_mass_in(i, 5) end do else ! 3D simulation do i = 1, num_fluids ! Loop through fluids write (i + 120, '(6X,6F24.12)') & nondim_time, & - c_mass(i, 1), & - c_mass(i, 2), & - c_mass(i, 3), & - c_mass(i, 4), & - c_mass(i, 5) + c_mass_in(i, 1), & + c_mass_in(i, 2), & + c_mass_in(i, 3), & + c_mass_in(i, 4), & + c_mass_in(i, 5) end do end if end if @@ -1071,7 +1070,7 @@ contains real(wp) :: max_pres real(wp), dimension(2) :: Re real(wp), dimension(6) :: tau_e - real(wp) :: G + real(wp) :: G_local real(wp) :: dyn_p, T real(wp) :: damage_state @@ -1152,7 +1151,7 @@ contains if (elasticity) then call s_convert_to_mixture_variables(q_cons_vf, j - 2, k, l, & rho, gamma, pi_inf, qv, & - Re, G, fluid_pp(:)%G) + Re, G_local, fluid_pp(:)%G) else call s_convert_to_mixture_variables(q_cons_vf, j - 2, k, l, & rho, gamma, pi_inf, qv) @@ -1166,7 +1165,7 @@ contains if (elasticity) then if (cont_damage) then damage_state = q_cons_vf(damage_idx)%sf(j - 2, k, l) - G = G*max((1._wp - damage_state), 0._wp) + G_local = G_local*max((1._wp - damage_state), 0._wp) end if call s_compute_pressure( & @@ -1174,7 +1173,7 @@ contains q_cons_vf(alf_idx)%sf(j - 2, k, l), & dyn_p, pi_inf, gamma, rho, qv, rhoYks(:), pres, T, & q_cons_vf(stress_idx%beg)%sf(j - 2, k, l), & - q_cons_vf(mom_idx%beg)%sf(j - 2, k, l), G) + q_cons_vf(mom_idx%beg)%sf(j - 2, k, l), G_local) else call s_compute_pressure( & q_cons_vf(1)%sf(j - 2, k, l), & @@ -1267,7 +1266,7 @@ contains ! Computing/Sharing necessary state variables call s_convert_to_mixture_variables(q_cons_vf, j - 2, k - 2, l, & rho, gamma, pi_inf, qv, & - Re, G, fluid_pp(:)%G) + Re, G_local, fluid_pp(:)%G) do s = 1, num_vels vel(s) = q_cons_vf(cont_idx%end + s)%sf(j - 2, k - 2, l)/rho end do @@ -1277,7 +1276,7 @@ contains if (elasticity) then if (cont_damage) then damage_state = q_cons_vf(damage_idx)%sf(j - 2, k - 2, l) - G = G*max((1._wp - damage_state), 0._wp) + G_local = G_local*max((1._wp - damage_state), 0._wp) end if call s_compute_pressure( & @@ -1288,7 +1287,7 @@ contains pres, & T, & q_cons_vf(stress_idx%beg)%sf(j - 2, k - 2, l), & - q_cons_vf(mom_idx%beg)%sf(j - 2, k - 2, l), G) + q_cons_vf(mom_idx%beg)%sf(j - 2, k - 2, l), G_local) else call s_compute_pressure(q_cons_vf(E_idx)%sf(j - 2, k - 2, l), & q_cons_vf(alf_idx)%sf(j - 2, k - 2, l), & @@ -1357,7 +1356,7 @@ contains ! Computing/Sharing necessary state variables call s_convert_to_mixture_variables(q_cons_vf, j - 2, k - 2, l - 2, & rho, gamma, pi_inf, qv, & - Re, G, fluid_pp(:)%G) + Re, G_local, fluid_pp(:)%G) do s = 1, num_vels vel(s) = q_cons_vf(cont_idx%end + s)%sf(j - 2, k - 2, l - 2)/rho end do @@ -1373,7 +1372,7 @@ contains if (elasticity) then if (cont_damage) then damage_state = q_cons_vf(damage_idx)%sf(j - 2, k - 2, l - 2) - G = G*max((1._wp - damage_state), 0._wp) + G_local = G_local*max((1._wp - damage_state), 0._wp) end if call s_compute_pressure( & @@ -1382,7 +1381,7 @@ contains dyn_p, pi_inf, gamma, rho, qv, & rhoYks, pres, T, & q_cons_vf(stress_idx%beg)%sf(j - 2, k - 2, l - 2), & - q_cons_vf(mom_idx%beg)%sf(j - 2, k - 2, l - 2), G) + q_cons_vf(mom_idx%beg)%sf(j - 2, k - 2, l - 2), G_local) else call s_compute_pressure(q_cons_vf(E_idx)%sf(j - 2, k - 2, l - 2), & q_cons_vf(alf_idx)%sf(j - 2, k - 2, l - 2), & diff --git a/src/simulation/m_fftw.fpp b/src/simulation/m_fftw.fpp index 87f612a76b..852fb90290 100644 --- a/src/simulation/m_fftw.fpp +++ b/src/simulation/m_fftw.fpp @@ -58,7 +58,6 @@ module m_fftw #else type(c_ptr) :: fwd_plan_gpu, bwd_plan_gpu #endif - integer :: ierr integer, allocatable :: gpu_fft_size(:), iembed(:), oembed(:) @@ -72,16 +71,18 @@ contains !! applying the Fourier filter in the azimuthal direction. impure subroutine s_initialize_fftw_module + integer :: ierr !< Generic flag used to identify and report GPU errors + ! Size of input array going into DFT real_size = p + 1 ! Size of output array coming out of DFT cmplx_size = (p + 1)/2 + 1 x_size = m + 1 - batch_size = x_size*sys_size #if defined(MFC_OpenACC) + rank = 1; istride = 1; ostride = 1 allocate (gpu_fft_size(1:rank), iembed(1:rank), oembed(1:rank)) @@ -134,11 +135,11 @@ contains real(c_double), pointer :: p_real(:) complex(c_double_complex), pointer :: p_cmplx(:), p_fltr_cmplx(:) integer :: i, j, k, l !< Generic loop iterators + integer :: ierr !< Generic flag used to identify and report GPU errors ! Restrict filter to processors that have cells adjacent to axis if (bc_y%beg >= 0) return #if defined(MFC_OpenACC) - $:GPU_PARALLEL_LOOP(collapse=3) do k = 1, sys_size do j = 0, m @@ -302,6 +303,7 @@ contains impure subroutine s_finalize_fftw_module #if defined(MFC_OpenACC) + integer :: ierr !< Generic flag used to identify and report GPU errors @:DEALLOCATE(data_real_gpu, data_fltr_cmplx_gpu, data_cmplx_gpu) #if defined(__PGI) diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index efe71526c6..69a7400062 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -244,8 +244,6 @@ module m_global_parameters integer :: mpi_info_int !> @} - integer, private :: ierr - !> @name Annotations of the structure of the state and flux vectors in terms of the !! size and the configuration of the system of equations to which they belong !> @{ @@ -1318,6 +1316,10 @@ contains !> Initializes parallel infrastructure impure subroutine s_initialize_parallel_io +#ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors +#endif + #:if not MFC_CASE_OPTIMIZATION num_dims = 1 + min(1, n) + min(1, p) diff --git a/src/simulation/m_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp index 628605a652..f4a24fba7a 100644 --- a/src/simulation/m_hyperelastic.fpp +++ b/src/simulation/m_hyperelastic.fpp @@ -103,11 +103,11 @@ contains real(wp), dimension(num_fluids) :: alpha_k, alpha_rho_k real(wp), dimension(2) :: Re real(wp) :: rho, gamma, pi_inf, qv - real(wp) :: G + real(wp) :: G_local integer :: j, k, l, i, r $:GPU_PARALLEL_LOOP(collapse=3, private='[alpha_K, alpha_rho_K, rho, & - & gamma, pi_inf, qv, G, Re, tensora, tensorb]') + & gamma, pi_inf, qv, G_local, Re, tensora, tensorb]') do l = 0, p do k = 0, n do j = 0, m @@ -118,12 +118,12 @@ contains end do ! If in simulation, use acc mixture subroutines call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha_k, & - alpha_rho_k, Re, G, Gs) + alpha_rho_k, Re, G_local, Gs) rho = max(rho, sgm_eps) - G = max(G, sgm_eps) - !if ( G <= verysmall ) G_K = 0._wp + G_local = max(G_local, sgm_eps) + !if ( G_local <= verysmall ) G_K = 0._wp - if (G > verysmall) then + if (G_local > verysmall) then $:GPU_LOOP(parallelism='[seq]') do i = 1, tensor_size tensora(i) = 0._wp @@ -190,13 +190,13 @@ contains btensor%vf(b_size)%sf(j, k, l) = tensorb(tensor_size) ! STEP 5a: updating the Cauchy stress primitive scalar field if (hyper_model == 1) then - call s_neoHookean_cauchy_solver(btensor%vf, q_prim_vf, G, j, k, l) + call s_neoHookean_cauchy_solver(btensor%vf, q_prim_vf, G_local, j, k, l) elseif (hyper_model == 2) then - call s_Mooney_Rivlin_cauchy_solver(btensor%vf, q_prim_vf, G, j, k, l) + call s_Mooney_Rivlin_cauchy_solver(btensor%vf, q_prim_vf, G_local, j, k, l) end if ! STEP 5b: updating the pressure field q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - & - G*q_prim_vf(xiend + 1)%sf(j, k, l)/gamma + G_local*q_prim_vf(xiend + 1)%sf(j, k, l)/gamma ! STEP 5c: updating the Cauchy stress conservative scalar field $:GPU_LOOP(parallelism='[seq]') do i = 1, b_size - 1 @@ -218,11 +218,11 @@ contains !! calculate the inverse of grad_xi to obtain F, F is a nxn tensor !! calculate the FFtranspose to obtain the btensor, btensor is nxn tensor !! btensor is symmetric, save the data space - pure subroutine s_neoHookean_cauchy_solver(btensor, q_prim_vf, G, j, k, l) + pure subroutine s_neoHookean_cauchy_solver(btensor_in, q_prim_vf, G_param, j, k, l) $:GPU_ROUTINE(parallelism='[seq]') type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf - type(scalar_field), dimension(b_size), intent(inout) :: btensor - real(wp), intent(in) :: G + type(scalar_field), dimension(b_size), intent(inout) :: btensor_in + real(wp), intent(in) :: G_param integer, intent(in) :: j, k, l real(wp) :: trace @@ -230,22 +230,22 @@ contains integer :: i !< Generic loop iterators ! tensor is the symmetric tensor & calculate the trace of the tensor - trace = btensor(1)%sf(j, k, l) + btensor(3)%sf(j, k, l) + btensor(6)%sf(j, k, l) + trace = btensor_in(1)%sf(j, k, l) + btensor_in(3)%sf(j, k, l) + btensor_in(6)%sf(j, k, l) ! calculate the deviatoric of the tensor #:for IJ in [1,3,6] - btensor(${IJ}$)%sf(j, k, l) = btensor(${IJ}$)%sf(j, k, l) - f13*trace + btensor_in(${IJ}$)%sf(j, k, l) = btensor_in(${IJ}$)%sf(j, k, l) - f13*trace #:endfor ! dividing by the jacobian for neo-Hookean model ! setting the tensor to the stresses for riemann solver $:GPU_LOOP(parallelism='[seq]') do i = 1, b_size - 1 q_prim_vf(strxb + i - 1)%sf(j, k, l) = & - G*btensor(i)%sf(j, k, l)/btensor(b_size)%sf(j, k, l) + G_param*btensor_in(i)%sf(j, k, l)/btensor_in(b_size)%sf(j, k, l) end do ! compute the invariant without the elastic modulus q_prim_vf(xiend + 1)%sf(j, k, l) = & - 0.5_wp*(trace - 3.0_wp)/btensor(b_size)%sf(j, k, l) + 0.5_wp*(trace - 3.0_wp)/btensor_in(b_size)%sf(j, k, l) end subroutine s_neoHookean_cauchy_solver @@ -257,11 +257,11 @@ contains !! calculate the inverse of grad_xi to obtain F, F is a nxn tensor !! calculate the FFtranspose to obtain the btensor, btensor is nxn tensor !! btensor is symmetric, save the data space - pure subroutine s_Mooney_Rivlin_cauchy_solver(btensor, q_prim_vf, G, j, k, l) + pure subroutine s_Mooney_Rivlin_cauchy_solver(btensor_in, q_prim_vf, G_param, j, k, l) $:GPU_ROUTINE(parallelism='[seq]') type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf - type(scalar_field), dimension(b_size), intent(inout) :: btensor - real(wp), intent(in) :: G + type(scalar_field), dimension(b_size), intent(inout) :: btensor_in + real(wp), intent(in) :: G_param integer, intent(in) :: j, k, l real(wp) :: trace @@ -270,23 +270,23 @@ contains !TODO Make this 1D and 2D capable ! tensor is the symmetric tensor & calculate the trace of the tensor - trace = btensor(1)%sf(j, k, l) + btensor(3)%sf(j, k, l) + btensor(6)%sf(j, k, l) + trace = btensor_in(1)%sf(j, k, l) + btensor_in(3)%sf(j, k, l) + btensor_in(6)%sf(j, k, l) ! calculate the deviatoric of the tensor - btensor(1)%sf(j, k, l) = btensor(1)%sf(j, k, l) - f13*trace - btensor(3)%sf(j, k, l) = btensor(3)%sf(j, k, l) - f13*trace - btensor(6)%sf(j, k, l) = btensor(6)%sf(j, k, l) - f13*trace + btensor_in(1)%sf(j, k, l) = btensor_in(1)%sf(j, k, l) - f13*trace + btensor_in(3)%sf(j, k, l) = btensor_in(3)%sf(j, k, l) - f13*trace + btensor_in(6)%sf(j, k, l) = btensor_in(6)%sf(j, k, l) - f13*trace ! dividing by the jacobian for neo-Hookean model ! setting the tensor to the stresses for riemann solver $:GPU_LOOP(parallelism='[seq]') do i = 1, b_size - 1 q_prim_vf(strxb + i - 1)%sf(j, k, l) = & - G*btensor(i)%sf(j, k, l)/btensor(b_size)%sf(j, k, l) + G_param*btensor_in(i)%sf(j, k, l)/btensor_in(b_size)%sf(j, k, l) end do ! compute the invariant without the elastic modulus q_prim_vf(xiend + 1)%sf(j, k, l) = & - 0.5_wp*(trace - 3.0_wp)/btensor(b_size)%sf(j, k, l) + 0.5_wp*(trace - 3.0_wp)/btensor_in(b_size)%sf(j, k, l) end subroutine s_Mooney_Rivlin_cauchy_solver diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index f9f12161b4..9b18f5b5fc 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -127,7 +127,7 @@ contains !! @param q_prim_vf Primitive variables !! @param pb Internal bubble pressure !! @param mv Mass of vapor in bubble - pure subroutine s_ibm_correct_state(q_cons_vf, q_prim_vf, pb, mv) + pure subroutine s_ibm_correct_state(q_cons_vf, q_prim_vf, pb_in, mv_in) type(scalar_field), & dimension(sys_size), & @@ -137,7 +137,7 @@ contains dimension(sys_size), & intent(INOUT) :: q_prim_vf !< Primitive Variables - real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), optional, intent(INOUT) :: pb, mv + real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), optional, intent(INOUT) :: pb_in, mv_in integer :: i, j, k, l, q, r!< Iterator variables integer :: patch_id !< Patch ID of ghost point @@ -198,7 +198,7 @@ contains else if (qbmm .and. .not. polytropic) then call s_interpolate_image_point(q_prim_vf, gp, & alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, & - r_IP, v_IP, pb_IP, mv_IP, nmom_IP, pb, mv, presb_IP, massv_IP) + r_IP, v_IP, pb_IP, mv_IP, nmom_IP, pb_in, mv_in, presb_IP, massv_IP) else call s_interpolate_image_point(q_prim_vf, gp, & alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP) @@ -297,8 +297,8 @@ contains if (.not. polytropic) then do q = 1, nb do r = 1, nnode - pb(j, k, l, r, q) = presb_IP((q - 1)*nnode + r) - mv(j, k, l, r, q) = massv_IP((q - 1)*nnode + r) + pb_in(j, k, l, r, q) = presb_IP((q - 1)*nnode + r) + mv_in(j, k, l, r, q) = massv_IP((q - 1)*nnode + r) end do end do end if @@ -336,11 +336,11 @@ contains !! @param ghost_points Ghost Points !! @param levelset Closest distance from each grid cell to IB !! @param levelset_norm Vector pointing in the direction of the closest distance - impure subroutine s_compute_image_points(ghost_points, levelset, levelset_norm) + impure subroutine s_compute_image_points(ghost_points_in, levelset_in, levelset_norm_in) - type(ghost_point), dimension(num_gps), intent(INOUT) :: ghost_points - type(levelset_field), intent(IN) :: levelset - type(levelset_norm_field), intent(IN) :: levelset_norm + type(ghost_point), dimension(num_gps), intent(INOUT) :: ghost_points_in + type(levelset_field), intent(IN) :: levelset_in + type(levelset_norm_field), intent(IN) :: levelset_norm_in real(wp) :: dist real(wp), dimension(3) :: norm @@ -357,7 +357,7 @@ contains integer :: index do q = 1, num_gps - gp = ghost_points(q) + gp = ghost_points_in(q) i = gp%loc(1) j = gp%loc(2) k = gp%loc(3) @@ -371,9 +371,9 @@ contains ! Calculate and store the precise location of the image point patch_id = gp%ib_patch_id - dist = abs(levelset%sf(i, j, k, patch_id)) - norm(:) = levelset_norm%sf(i, j, k, patch_id, :) - ghost_points(q)%ip_loc(:) = physical_loc(:) + 2*dist*norm(:) + dist = abs(levelset_in%sf(i, j, k, patch_id)) + norm(:) = levelset_norm_in%sf(i, j, k, patch_id, :) + ghost_points_in(q)%ip_loc(:) = physical_loc(:) + 2*dist*norm(:) ! Find the closest grid point to the image point do dim = 1, num_dims @@ -391,7 +391,7 @@ contains end if if (f_approx_equal(norm(dim), 0._wp)) then - ghost_points(q)%ip_grid(dim) = ghost_points(q)%loc(dim) + ghost_points_in(q)%ip_grid(dim) = ghost_points_in(q)%loc(dim) else if (norm(dim) > 0) then dir = 1 @@ -399,18 +399,18 @@ contains dir = -1 end if - index = ghost_points(q)%loc(dim) - temp_loc = ghost_points(q)%ip_loc(dim) + index = ghost_points_in(q)%loc(dim) + temp_loc = ghost_points_in(q)%ip_loc(dim) do while ((temp_loc < s_cc(index) & .or. temp_loc > s_cc(index + 1)) & .and. (index >= 0 .and. index <= bound)) index = index + dir end do - ghost_points(q)%ip_grid(dim) = index - if (ghost_points(q)%DB(dim) == -1) then - ghost_points(q)%ip_grid(dim) = ghost_points(q)%loc(dim) + 1 - else if (ghost_points(q)%DB(dim) == 1) then - ghost_points(q)%ip_grid(dim) = ghost_points(q)%loc(dim) - 1 + ghost_points_in(q)%ip_grid(dim) = index + if (ghost_points_in(q)%DB(dim) == -1) then + ghost_points_in(q)%ip_grid(dim) = ghost_points_in(q)%loc(dim) + 1 + else if (ghost_points_in(q)%DB(dim) == 1) then + ghost_points_in(q)%ip_grid(dim) = ghost_points_in(q)%loc(dim) - 1 end if end if end do @@ -420,10 +420,10 @@ contains !> Function that finds the number of ghost points, used for allocating !! memory. - pure subroutine s_find_num_ghost_points(num_gps, num_inner_gps) + pure subroutine s_find_num_ghost_points(num_gps_out, num_inner_gps_out) - integer, intent(out) :: num_gps - integer, intent(out) :: num_inner_gps + integer, intent(out) :: num_gps_out + integer, intent(out) :: num_inner_gps_out integer, dimension(2*gp_layers + 1, 2*gp_layers + 1) & :: subsection_2D @@ -431,8 +431,8 @@ contains :: subsection_3D integer :: i, j, k!< Iterator variables - num_gps = 0 - num_inner_gps = 0 + num_gps_out = 0 + num_inner_gps_out = 0 do i = 0, m do j = 0, n @@ -442,9 +442,9 @@ contains i - gp_layers:i + gp_layers, & j - gp_layers:j + gp_layers, 0) if (any(subsection_2D == 0)) then - num_gps = num_gps + 1 + num_gps_out = num_gps_out + 1 else - num_inner_gps = num_inner_gps + 1 + num_inner_gps_out = num_inner_gps_out + 1 end if end if else @@ -455,9 +455,9 @@ contains j - gp_layers:j + gp_layers, & k - gp_layers:k + gp_layers) if (any(subsection_3D == 0)) then - num_gps = num_gps + 1 + num_gps_out = num_gps_out + 1 else - num_inner_gps = num_inner_gps + 1 + num_inner_gps_out = num_inner_gps_out + 1 end if end if end do @@ -468,10 +468,10 @@ contains end subroutine s_find_num_ghost_points !> Function that finds the ghost points - pure subroutine s_find_ghost_points(ghost_points, inner_points) + pure subroutine s_find_ghost_points(ghost_points_in, inner_points_in) - type(ghost_point), dimension(num_gps), intent(INOUT) :: ghost_points - type(ghost_point), dimension(num_inner_gps), intent(INOUT) :: inner_points + type(ghost_point), dimension(num_gps), intent(INOUT) :: ghost_points_in + type(ghost_point), dimension(num_inner_gps), intent(INOUT) :: inner_points_in integer, dimension(2*gp_layers + 1, 2*gp_layers + 1) & :: subsection_2D integer, dimension(2*gp_layers + 1, 2*gp_layers + 1, 2*gp_layers + 1) & @@ -491,37 +491,37 @@ contains i - gp_layers:i + gp_layers, & j - gp_layers:j + gp_layers, 0) if (any(subsection_2D == 0)) then - ghost_points(count)%loc = [i, j, 0] + ghost_points_in(count)%loc = [i, j, 0] patch_id = ib_markers%sf(i, j, 0) - ghost_points(count)%ib_patch_id = & + ghost_points_in(count)%ib_patch_id = & patch_id - ghost_points(count)%slip = patch_ib(patch_id)%slip + ghost_points_in(count)%slip = patch_ib(patch_id)%slip ! ghost_points(count)%rank = proc_rank if ((x_cc(i) - dx(i)) < x_domain%beg) then - ghost_points(count)%DB(1) = -1 + ghost_points_in(count)%DB(1) = -1 else if ((x_cc(i) + dx(i)) > x_domain%end) then - ghost_points(count)%DB(1) = 1 + ghost_points_in(count)%DB(1) = 1 else - ghost_points(count)%DB(1) = 0 + ghost_points_in(count)%DB(1) = 0 end if if ((y_cc(j) - dy(j)) < y_domain%beg) then - ghost_points(count)%DB(2) = -1 + ghost_points_in(count)%DB(2) = -1 else if ((y_cc(j) + dy(j)) > y_domain%end) then - ghost_points(count)%DB(2) = 1 + ghost_points_in(count)%DB(2) = 1 else - ghost_points(count)%DB(2) = 0 + ghost_points_in(count)%DB(2) = 0 end if count = count + 1 else - inner_points(count_i)%loc = [i, j, 0] + inner_points_in(count_i)%loc = [i, j, 0] patch_id = ib_markers%sf(i, j, 0) - inner_points(count_i)%ib_patch_id = & + inner_points_in(count_i)%ib_patch_id = & patch_id - inner_points(count_i)%slip = patch_ib(patch_id)%slip + inner_points_in(count_i)%slip = patch_ib(patch_id)%slip count_i = count_i + 1 end if @@ -534,43 +534,43 @@ contains j - gp_layers:j + gp_layers, & k - gp_layers:k + gp_layers) if (any(subsection_3D == 0)) then - ghost_points(count)%loc = [i, j, k] + ghost_points_in(count)%loc = [i, j, k] patch_id = ib_markers%sf(i, j, k) - ghost_points(count)%ib_patch_id = & + ghost_points_in(count)%ib_patch_id = & ib_markers%sf(i, j, k) - ghost_points(count)%slip = patch_ib(patch_id)%slip + ghost_points_in(count)%slip = patch_ib(patch_id)%slip if ((x_cc(i) - dx(i)) < x_domain%beg) then - ghost_points(count)%DB(1) = -1 + ghost_points_in(count)%DB(1) = -1 else if ((x_cc(i) + dx(i)) > x_domain%end) then - ghost_points(count)%DB(1) = 1 + ghost_points_in(count)%DB(1) = 1 else - ghost_points(count)%DB(1) = 0 + ghost_points_in(count)%DB(1) = 0 end if if ((y_cc(j) - dy(j)) < y_domain%beg) then - ghost_points(count)%DB(2) = -1 + ghost_points_in(count)%DB(2) = -1 else if ((y_cc(j) + dy(j)) > y_domain%end) then - ghost_points(count)%DB(2) = 1 + ghost_points_in(count)%DB(2) = 1 else - ghost_points(count)%DB(2) = 0 + ghost_points_in(count)%DB(2) = 0 end if if ((z_cc(k) - dz(k)) < z_domain%beg) then - ghost_points(count)%DB(3) = -1 + ghost_points_in(count)%DB(3) = -1 else if ((z_cc(k) + dz(k)) > z_domain%end) then - ghost_points(count)%DB(3) = 1 + ghost_points_in(count)%DB(3) = 1 else - ghost_points(count)%DB(3) = 0 + ghost_points_in(count)%DB(3) = 0 end if count = count + 1 else - inner_points(count_i)%loc = [i, j, k] + inner_points_in(count_i)%loc = [i, j, k] patch_id = ib_markers%sf(i, j, k) - inner_points(count_i)%ib_patch_id = & + inner_points_in(count_i)%ib_patch_id = & ib_markers%sf(i, j, k) - inner_points(count_i)%slip = patch_ib(patch_id)%slip + inner_points_in(count_i)%slip = patch_ib(patch_id)%slip count_i = count_i + 1 end if @@ -583,9 +583,9 @@ contains end subroutine s_find_ghost_points !> Function that computes the interpolation coefficients of image points - pure subroutine s_compute_interpolation_coeffs(ghost_points) + pure subroutine s_compute_interpolation_coeffs(ghost_points_in) - type(ghost_point), dimension(num_gps), intent(INOUT) :: ghost_points + type(ghost_point), dimension(num_gps), intent(INOUT) :: ghost_points_in real(wp), dimension(2, 2, 2) :: dist real(wp), dimension(2, 2, 2) :: alpha @@ -600,7 +600,7 @@ contains ! 2D if (p <= 0) then do i = 1, num_gps - gp = ghost_points(i) + gp = ghost_points_in(i) ! Get the interpolation points i1 = gp%ip_grid(1); i2 = i1 + 1 j1 = gp%ip_grid(2); j2 = j1 + 1 @@ -647,12 +647,12 @@ contains end if end if - ghost_points(i)%interp_coeffs = interp_coeffs + ghost_points_in(i)%interp_coeffs = interp_coeffs end do else do i = 1, num_gps - gp = ghost_points(i) + gp = ghost_points_in(i) ! Get the interpolation points i1 = gp%ip_grid(1); i2 = i1 + 1 j1 = gp%ip_grid(2); j2 = j1 + 1 @@ -729,7 +729,7 @@ contains end if end if - ghost_points(i)%interp_coeffs = interp_coeffs + ghost_points_in(i)%interp_coeffs = interp_coeffs end do end if @@ -737,13 +737,13 @@ contains !> Function that uses the interpolation coefficients and the current state !! at the cell centers in order to estimate the state at the image point - pure subroutine s_interpolate_image_point(q_prim_vf, gp, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, r_IP, v_IP, pb_IP, mv_IP, nmom_IP, pb, mv, presb_IP, massv_IP) + pure subroutine s_interpolate_image_point(q_prim_vf, gp, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, r_IP, v_IP, pb_IP, mv_IP, nmom_IP, pb_in, mv_in, presb_IP, massv_IP) $:GPU_ROUTINE(parallelism='[seq]') type(scalar_field), & dimension(sys_size), & intent(IN) :: q_prim_vf !< Primitive Variables - real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(INOUT) :: pb, mv + real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(INOUT) :: pb_in, mv_in type(ghost_point), intent(IN) :: gp real(wp), intent(INOUT) :: pres_IP @@ -843,8 +843,8 @@ contains if (.not. polytropic) then do q = 1, nb do l = 1, nnode - presb_IP((q - 1)*nnode + l) = presb_IP((q - 1)*nnode + l) + coeff*pb(i, j, k, l, q) - massv_IP((q - 1)*nnode + l) = massv_IP((q - 1)*nnode + l) + coeff*mv(i, j, k, l, q) + presb_IP((q - 1)*nnode + l) = presb_IP((q - 1)*nnode + l) + coeff*pb_in(i, j, k, l, q) + massv_IP((q - 1)*nnode + l) = massv_IP((q - 1)*nnode + l) + coeff*mv_in(i, j, k, l, q) end do end do end if diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index 07a21e1c96..3054f6f66c 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -42,11 +42,6 @@ module m_mpi_proxy !! immersed boundary markers, for a single computational domain boundary !! at the time, from the relevant neighboring processor. - !> @name Generic flags used to identify and report MPI errors - !> @{ - integer, private :: ierr - !> @} - integer :: i_halo_size $:GPU_DECLARE(create='[i_halo_size]') @@ -88,6 +83,7 @@ contains #ifdef MFC_MPI integer :: i, j !< Generic loop iterator + integer :: ierr !< Generic flag used to identify and report MPI errors call MPI_BCAST(case_dir, len(case_dir), MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr) @@ -259,6 +255,7 @@ contains integer :: pack_offset, unpack_offset #ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors call nvtxStartRange("IB-MARKER-COMM-PACKBUF") @@ -396,6 +393,7 @@ contains real(wp), intent(inout), dimension(1:num_freq) :: phi_rn #ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors call MPI_BCAST(phi_rn, num_freq, mpi_p, 0, MPI_COMM_WORLD, ierr) #endif diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index b733cfa2dd..ae66f00011 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -614,15 +614,15 @@ contains end subroutine s_initialize_rhs_module - impure subroutine s_compute_rhs(q_cons_vf, q_T_sf, q_prim_vf, bc_type, rhs_vf, pb, rhs_pb, mv, rhs_mv, t_step, time_avg, stage) + impure subroutine s_compute_rhs(q_cons_vf, q_T_sf, q_prim_vf, bc_type, rhs_vf, pb_in, rhs_pb, mv_in, rhs_mv, t_step, time_avg, stage) type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf type(scalar_field), intent(inout) :: q_T_sf type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf type(integer_field), dimension(1:num_dims, -1:1), intent(in) :: bc_type type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf - real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, rhs_pb - real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: mv, rhs_mv + real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb_in, rhs_pb + real(wp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: mv_in, rhs_mv integer, intent(in) :: t_step real(wp), intent(inout) :: time_avg integer, intent(in) :: stage @@ -672,7 +672,7 @@ contains if (igr) then call nvtxStartRange("RHS-COMMUNICATION") - call s_populate_variables_buffers(bc_type, q_cons_vf, pb, mv) + call s_populate_variables_buffers(bc_type, q_cons_vf, pb_in, mv_in) call nvtxEndRange else call nvtxStartRange("RHS-CONVERT") @@ -684,7 +684,7 @@ contains call nvtxEndRange call nvtxStartRange("RHS-COMMUNICATION") - call s_populate_variables_buffers(bc_type, q_prim_qp%vf, pb, mv) + call s_populate_variables_buffers(bc_type, q_prim_qp%vf, pb_in, mv_in) call nvtxEndRange end if @@ -698,7 +698,7 @@ contains if (t_step == t_step_stop) return end if - if (qbmm) call s_mom_inv(q_cons_qp%vf, q_prim_qp%vf, mom_sp, mom_3d, pb, rhs_pb, mv, rhs_mv, idwbuff(1), idwbuff(2), idwbuff(3)) + if (qbmm) call s_mom_inv(q_cons_qp%vf, q_prim_qp%vf, mom_sp, mom_3d, pb_in, rhs_pb, mv_in, rhs_mv, idwbuff(1), idwbuff(2), idwbuff(3)) if (viscous .and. .not. igr) then call nvtxStartRange("RHS-VISCOUS") @@ -716,7 +716,7 @@ contains if (surface_tension) then call nvtxStartRange("RHS-SURFACE-TENSION") - call s_get_capilary(q_prim_qp%vf, bc_type) + call s_get_capillary(q_prim_qp%vf, bc_type) call nvtxEndRange end if @@ -892,7 +892,7 @@ contains q_prim_qp%vf, & rhs_vf, & flux_n(id)%vf, & - pb, & + pb_in, & rhs_pb) call nvtxEndRange end if @@ -939,7 +939,8 @@ contains call s_compute_bubble_EE_source( & q_cons_qp%vf(1:sys_size), & q_prim_qp%vf(1:sys_size), & - rhs_vf) + rhs_vf, & + divu) call nvtxEndRange end if @@ -1468,13 +1469,13 @@ contains end subroutine s_compute_advection_source_term - subroutine s_compute_additional_physics_rhs(idir, q_prim_vf, rhs_vf, flux_src_n, & + subroutine s_compute_additional_physics_rhs(idir, q_prim_vf, rhs_vf, flux_src_n_in, & dq_prim_dx_vf, dq_prim_dy_vf, dq_prim_dz_vf) integer, intent(in) :: idir type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf - type(scalar_field), dimension(sys_size), intent(in) :: flux_src_n + type(scalar_field), dimension(sys_size), intent(in) :: flux_src_n_in type(scalar_field), dimension(sys_size), intent(in) :: dq_prim_dx_vf, dq_prim_dy_vf, dq_prim_dz_vf integer :: i, j, k, l @@ -1489,8 +1490,8 @@ contains rhs_vf(c_idx)%sf(j, k, l) = & rhs_vf(c_idx)%sf(j, k, l) + 1._wp/dx(j)* & q_prim_vf(c_idx)%sf(j, k, l)* & - (flux_src_n(advxb)%sf(j, k, l) - & - flux_src_n(advxb)%sf(j - 1, k, l)) + (flux_src_n_in(advxb)%sf(j, k, l) - & + flux_src_n_in(advxb)%sf(j - 1, k, l)) end do end do end do @@ -1504,8 +1505,8 @@ contains do i = momxb, E_idx rhs_vf(i)%sf(j, k, l) = & rhs_vf(i)%sf(j, k, l) + 1._wp/dx(j)* & - (flux_src_n(i)%sf(j - 1, k, l) & - - flux_src_n(i)%sf(j, k, l)) + (flux_src_n_in(i)%sf(j - 1, k, l) & + - flux_src_n_in(i)%sf(j, k, l)) end do end do end do @@ -1521,8 +1522,8 @@ contains rhs_vf(c_idx)%sf(j, k, l) = & rhs_vf(c_idx)%sf(j, k, l) + 1._wp/dy(k)* & q_prim_vf(c_idx)%sf(j, k, l)* & - (flux_src_n(advxb)%sf(j, k, l) - & - flux_src_n(advxb)%sf(j, k - 1, l)) + (flux_src_n_in(advxb)%sf(j, k, l) - & + flux_src_n_in(advxb)%sf(j, k - 1, l)) end do end do end do @@ -1569,8 +1570,8 @@ contains do i = momxb, E_idx rhs_vf(i)%sf(j, k, l) = & rhs_vf(i)%sf(j, k, l) + 1._wp/dy(k)* & - (flux_src_n(i)%sf(j, k - 1, l) & - - flux_src_n(i)%sf(j, k, l)) + (flux_src_n_in(i)%sf(j, k - 1, l) & + - flux_src_n_in(i)%sf(j, k, l)) end do end do end do @@ -1585,8 +1586,8 @@ contains do i = momxb, E_idx rhs_vf(i)%sf(j, k, l) = & rhs_vf(i)%sf(j, k, l) + 1._wp/dy(k)* & - (flux_src_n(i)%sf(j, k - 1, l) & - - flux_src_n(i)%sf(j, k, l)) + (flux_src_n_in(i)%sf(j, k - 1, l) & + - flux_src_n_in(i)%sf(j, k, l)) end do end do end do @@ -1606,8 +1607,8 @@ contains do i = momxb, E_idx rhs_vf(i)%sf(j, k, l) = & rhs_vf(i)%sf(j, k, l) - 5.e-1_wp/y_cc(k)* & - (flux_src_n(i)%sf(j, k - 1, l) & - + flux_src_n(i)%sf(j, k, l)) + (flux_src_n_in(i)%sf(j, k - 1, l) & + + flux_src_n_in(i)%sf(j, k, l)) end do end do end do @@ -1636,8 +1637,8 @@ contains do i = momxb, E_idx rhs_vf(i)%sf(j, k, l) = & rhs_vf(i)%sf(j, k, l) - 5.e-1_wp/y_cc(k)* & - (flux_src_n(i)%sf(j, k - 1, l) & - + flux_src_n(i)%sf(j, k, l)) + (flux_src_n_in(i)%sf(j, k - 1, l) & + + flux_src_n_in(i)%sf(j, k, l)) end do end do end do @@ -1656,8 +1657,8 @@ contains rhs_vf(c_idx)%sf(j, k, l) = & rhs_vf(c_idx)%sf(j, k, l) + 1._wp/dz(l)* & q_prim_vf(c_idx)%sf(j, k, l)* & - (flux_src_n(advxb)%sf(j, k, l) - & - flux_src_n(advxb)%sf(j, k, l - 1)) + (flux_src_n_in(advxb)%sf(j, k, l) - & + flux_src_n_in(advxb)%sf(j, k, l - 1)) end do end do end do @@ -1671,8 +1672,8 @@ contains do i = momxb, E_idx rhs_vf(i)%sf(j, k, l) = & rhs_vf(i)%sf(j, k, l) + 1._wp/dz(l)* & - (flux_src_n(i)%sf(j, k, l - 1) & - - flux_src_n(i)%sf(j, k, l)) + (flux_src_n_in(i)%sf(j, k, l - 1) & + - flux_src_n_in(i)%sf(j, k, l)) end do end do end do @@ -1685,13 +1686,13 @@ contains do j = 0, m rhs_vf(momxb + 1)%sf(j, k, l) = & rhs_vf(momxb + 1)%sf(j, k, l) + 5.e-1_wp* & - (flux_src_n(momxe)%sf(j, k, l - 1) & - + flux_src_n(momxe)%sf(j, k, l)) + (flux_src_n_in(momxe)%sf(j, k, l - 1) & + + flux_src_n_in(momxe)%sf(j, k, l)) rhs_vf(momxe)%sf(j, k, l) = & rhs_vf(momxe)%sf(j, k, l) - 5.e-1_wp* & - (flux_src_n(momxb + 1)%sf(j, k, l - 1) & - + flux_src_n(momxb + 1)%sf(j, k, l)) + (flux_src_n_in(momxb + 1)%sf(j, k, l - 1) & + + flux_src_n_in(momxb + 1)%sf(j, k, l)) end do end do end do @@ -1748,7 +1749,6 @@ contains if (n > 0) then if (p > 0) then - call s_weno(v_vf(iv%beg:iv%end), & vL_x(:, :, :, iv%beg:iv%end), vL_y(:, :, :, iv%beg:iv%end), vL_z(:, :, :, iv%beg:iv%end), vR_x(:, :, :, iv%beg:iv%end), vR_y(:, :, :, iv%beg:iv%end), vR_z(:, :, :, iv%beg:iv%end), & weno_dir, & @@ -1760,7 +1760,6 @@ contains is1, is2, is3) end if else - call s_weno(v_vf(iv%beg:iv%end), & vL_x(:, :, :, iv%beg:iv%end), vL_y(:, :, :, :), vL_z(:, :, :, :), vR_x(:, :, :, iv%beg:iv%end), vR_y(:, :, :, :), vR_z(:, :, :, :), & weno_dir, & diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 1b77e17bfd..8196403564 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -37,7 +37,7 @@ module m_riemann_solvers use m_bubbles_EE - use m_surface_tension !< To get the capilary fluxes + use m_surface_tension !< To get the capillary fluxes use m_helper_basic !< Functions to compare floating point numbers @@ -2887,7 +2887,7 @@ contains end if if (surface_tension) then - call s_compute_capilary_source_flux( & + call s_compute_capillary_source_flux( & vel_src_rsx_vf, & vel_src_rsy_vf, & vel_src_rsz_vf, & diff --git a/src/simulation/m_sim_helpers.fpp b/src/simulation/m_sim_helpers.fpp index 2105a43d0b..565524e80b 100644 --- a/src/simulation/m_sim_helpers.fpp +++ b/src/simulation/m_sim_helpers.fpp @@ -101,7 +101,7 @@ contains real(wp), dimension(2), intent(inout) :: Re real(wp), dimension(num_fluids) :: alpha_rho, Gs - real(wp) :: qv, E, G + real(wp) :: qv, E, G_local integer :: i @@ -129,7 +129,7 @@ contains if (elasticity) then call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, & - alpha_rho, Re, G, Gs) + alpha_rho, Re, G_local, Gs) elseif (bubbles_euler) then call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re) else @@ -164,7 +164,7 @@ contains ! Adjust energy for hyperelasticity if (hyperelasticity) then - E = E + G*q_prim_vf(xiend + 1)%sf(j, k, l) + E = E + G_local*q_prim_vf(xiend + 1)%sf(j, k, l) end if H = (E + pres)/rho diff --git a/src/simulation/m_surface_tension.fpp b/src/simulation/m_surface_tension.fpp index 67c31cf74e..5d23d8e4c3 100644 --- a/src/simulation/m_surface_tension.fpp +++ b/src/simulation/m_surface_tension.fpp @@ -21,8 +21,8 @@ module m_surface_tension implicit none private; public :: s_initialize_surface_tension_module, & - s_compute_capilary_source_flux, & - s_get_capilary, & + s_compute_capillary_source_flux, & + s_get_capillary, & s_finalize_surface_tension_module !> @name color function gradient components and magnitude @@ -65,7 +65,7 @@ contains end if end subroutine s_initialize_surface_tension_module - pure subroutine s_compute_capilary_source_flux( & + pure subroutine s_compute_capillary_source_flux( & vSrc_rsx_vf, vSrc_rsy_vf, vSrc_rsz_vf, & flux_src_vf, & id, isx, isy, isz) @@ -111,7 +111,7 @@ contains normW = (normWL + normWR)/2._wp if (normW > capillary_cutoff) then - @:compute_capilary_stress_tensor() + @:compute_capillary_stress_tensor() do i = 1, num_dims @@ -158,7 +158,7 @@ contains normW = (normWL + normWR)/2._wp if (normW > capillary_cutoff) then - @:compute_capilary_stress_tensor() + @:compute_capillary_stress_tensor() do i = 1, num_dims @@ -205,7 +205,7 @@ contains normW = (normWL + normWR)/2._wp if (normW > capillary_cutoff) then - @:compute_capilary_stress_tensor() + @:compute_capillary_stress_tensor() do i = 1, num_dims @@ -226,9 +226,9 @@ contains end if - end subroutine s_compute_capilary_source_flux + end subroutine s_compute_capillary_source_flux - impure subroutine s_get_capilary(q_prim_vf, bc_type) + impure subroutine s_get_capillary(q_prim_vf, bc_type) type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf type(integer_field), dimension(1:num_dims, -1:1), intent(in) :: bc_type @@ -301,7 +301,7 @@ contains call s_reconstruct_cell_boundary_values_capillary(c_divs, gL_x, gL_y, gL_z, gR_x, gR_y, gR_z, i) end do - end subroutine s_get_capilary + end subroutine s_get_capillary subroutine s_reconstruct_cell_boundary_values_capillary(v_vf, vL_x, vL_y, vL_z, vR_x, vR_y, vR_z, & norm_dir) diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 725fa15a7d..bd9bd1ef54 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -955,7 +955,7 @@ contains if (bubbles_euler) then - call s_compute_bubble_EE_source(q_cons_ts(1)%vf, q_prim_vf, rhs_vf) + call s_compute_bubble_EE_source(q_cons_ts(1)%vf, q_prim_vf, rhs_vf, divu) call s_comp_alpha_from_n(q_cons_ts(1)%vf) elseif (bubbles_lagrange) then @@ -1036,18 +1036,18 @@ contains !> This subroutine applies the body forces source term at each !! Runge-Kutta stage - subroutine s_apply_bodyforces(q_cons_vf, q_prim_vf, rhs_vf, ldt) + subroutine s_apply_bodyforces(q_cons_vf, q_prim_vf_in, rhs_vf_in, ldt) type(scalar_field), dimension(1:sys_size), intent(inout) :: q_cons_vf - type(scalar_field), dimension(1:sys_size), intent(in) :: q_prim_vf - type(scalar_field), dimension(1:sys_size), intent(inout) :: rhs_vf + type(scalar_field), dimension(1:sys_size), intent(in) :: q_prim_vf_in + type(scalar_field), dimension(1:sys_size), intent(inout) :: rhs_vf_in real(wp), intent(in) :: ldt !< local dt integer :: i, j, k, l call nvtxStartRange("RHS-BODYFORCES") - call s_compute_body_forces_rhs(q_prim_vf, q_cons_vf, rhs_vf) + call s_compute_body_forces_rhs(q_prim_vf_in, q_cons_vf, rhs_vf_in) $:GPU_PARALLEL_LOOP(collapse=4) do i = momxb, E_idx @@ -1055,7 +1055,7 @@ contains do k = 0, n do j = 0, m q_cons_vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l) + & - ldt*rhs_vf(i)%sf(j, k, l) + ldt*rhs_vf_in(i)%sf(j, k, l) end do end do end do