From d57312425676c5525ac5a2e8ee7b04bb3f98afbf Mon Sep 17 00:00:00 2001 From: AppleChen17 Date: Sat, 17 Oct 2026 01:35:45 +0800 Subject: [PATCH] Collapse only spatial loops in per-variable copies Eight loops copy or difference every variable of a field using collapse=4 over the variable and three spatial dimensions. Move the variable loop innermost and execute it sequentially, collapsing only the three spatial loops. In the flux-difference loops, the inverse grid spacing is also computed outside the variable loop. Affected routines: - m_rhs: s_compute_rhs, s_compute_advection_source_term - m_riemann_state: s_finalize_riemann_solver, s_finalize_riemann_solver_hatR - m_weno: s_pack_weno_input_arr The per-variable numerical expressions and spatial loop order are preserved; outputs are bit-identical to the baseline in our tests. --- src/simulation/m_rhs.fpp | 50 ++++++++++++++++-------------- src/simulation/m_riemann_state.fpp | 33 +++++++++++--------- src/simulation/m_weno.fpp | 11 ++++--- 3 files changed, 51 insertions(+), 43 deletions(-) diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index 7e74b33ed2..66c0fe934b 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -599,11 +599,12 @@ contains if (.not. igr) then ! Association/Population of Working Variables - $:GPU_PARALLEL_LOOP(private='[i, j, k, l]', collapse=4) - do i = 1, sys_size - do l = idwbuff(3)%beg, idwbuff(3)%end - do k = idwbuff(2)%beg, idwbuff(2)%end - do j = idwbuff(1)%beg, idwbuff(1)%end + $:GPU_PARALLEL_LOOP(private='[i, j, k, l]', collapse=3) + do l = idwbuff(3)%beg, idwbuff(3)%end + do k = idwbuff(2)%beg, idwbuff(2)%end + do j = idwbuff(1)%beg, idwbuff(1)%end + $:GPU_LOOP(parallelism='[seq]') + do i = 1, sys_size q_cons_qp%vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l) end do end do @@ -1124,12 +1125,13 @@ contains end if if (hypo_nc_mode /= hypo_nc_mode_dual_pass) then - $:GPU_PARALLEL_LOOP(collapse=4,private='[j, k_loop, l_loop, q_loop, inv_ds, flux_face1, flux_face2]') - do j = 1, sys_size - do q_loop = 0, p - do l_loop = 0, n - do k_loop = 0, m - inv_ds = 1._wp/dx(k_loop) + $:GPU_PARALLEL_LOOP(collapse=3,private='[j, k_loop, l_loop, q_loop, inv_ds, flux_face1, flux_face2]') + do q_loop = 0, p + do l_loop = 0, n + do k_loop = 0, m + inv_ds = 1._wp/dx(k_loop) + $:GPU_LOOP(parallelism='[seq]') + do j = 1, sys_size flux_face1 = flux_n(1)%vf(j)%sf(k_loop - 1, l_loop, q_loop) flux_face2 = flux_n(1)%vf(j)%sf(k_loop, l_loop, q_loop) rhs_vf(j)%sf(k_loop, l_loop, q_loop) = inv_ds*(flux_face1 - flux_face2) @@ -1200,12 +1202,13 @@ contains end if if (hypo_nc_mode /= hypo_nc_mode_dual_pass) then - $:GPU_PARALLEL_LOOP(collapse=4,private='[j, k, l, q, inv_ds, flux_face1, flux_face2]') - do j = 1, sys_size - do l = 0, p - do k = 0, n - do q = 0, m - inv_ds = 1._wp/dy(k) + $:GPU_PARALLEL_LOOP(collapse=3,private='[j, k, l, q, inv_ds, flux_face1, flux_face2]') + do l = 0, p + do k = 0, n + do q = 0, m + inv_ds = 1._wp/dy(k) + $:GPU_LOOP(parallelism='[seq]') + do j = 1, sys_size flux_face1 = flux_n(2)%vf(j)%sf(q, k - 1, l) flux_face2 = flux_n(2)%vf(j)%sf(q, k, l) rhs_vf(j)%sf(q, k, l) = rhs_vf(j)%sf(q, k, l) + inv_ds*(flux_face1 - flux_face2) @@ -1362,12 +1365,13 @@ contains $:END_GPU_PARALLEL_LOOP() else ! Cartesian Coordinates if (hypo_nc_mode /= hypo_nc_mode_dual_pass) then - $:GPU_PARALLEL_LOOP(collapse=4,private='[j, k, l, q, inv_ds, flux_face1, flux_face2]') - do j = 1, sys_size - do k = 0, p - do q = 0, n - do l = 0, m - inv_ds = 1._wp/dz(k) + $:GPU_PARALLEL_LOOP(collapse=3,private='[j, k, l, q, inv_ds, flux_face1, flux_face2]') + do k = 0, p + do q = 0, n + do l = 0, m + inv_ds = 1._wp/dz(k) + $:GPU_LOOP(parallelism='[seq]') + do j = 1, sys_size flux_face1 = flux_n(3)%vf(j)%sf(l, q, k - 1) flux_face2 = flux_n(3)%vf(j)%sf(l, q, k) rhs_vf(j)%sf(l, q, k) = rhs_vf(j)%sf(l, q, k) + inv_ds*(flux_face1 - flux_face2) diff --git a/src/simulation/m_riemann_state.fpp b/src/simulation/m_riemann_state.fpp index aeb34c8a9f..1c589d4f7f 100644 --- a/src/simulation/m_riemann_state.fpp +++ b/src/simulation/m_riemann_state.fpp @@ -1345,11 +1345,12 @@ contains ! Reshaping Outputted Data in y-direction if (norm_dir == 2) then - $:GPU_PARALLEL_LOOP(collapse=4) - do i = 1, sys_size - do l = is3%beg, is3%end - do j = is1%beg, is1%end - do k = is2%beg, is2%end + $:GPU_PARALLEL_LOOP(collapse=3) + do l = is3%beg, is3%end + do j = is1%beg, is1%end + do k = is2%beg, is2%end + $:GPU_LOOP(parallelism='[seq]') + do i = 1, sys_size flux_vf(i)%sf(k, j, l) = flux_${INFIX}$rsx_vf(k, j, l, i) end do end do @@ -1401,11 +1402,12 @@ contains #:endif ! Reshaping Outputted Data in z-direction else if (norm_dir == 3) then - $:GPU_PARALLEL_LOOP(collapse=4) - do i = 1, sys_size - do j = is1%beg, is1%end - do k = is2%beg, is2%end - do l = is3%beg, is3%end + $:GPU_PARALLEL_LOOP(collapse=3) + do j = is1%beg, is1%end + do k = is2%beg, is2%end + do l = is3%beg, is3%end + $:GPU_LOOP(parallelism='[seq]') + do i = 1, sys_size flux_vf(i)%sf(l, k, j) = flux_${INFIX}$rsx_vf(l, k, j, i) end do end do @@ -1455,11 +1457,12 @@ contains end if #:endif else if (norm_dir == 1) then - $:GPU_PARALLEL_LOOP(collapse=4) - do i = 1, sys_size - do l = is3%beg, is3%end - do k = is2%beg, is2%end - do j = is1%beg, is1%end + $:GPU_PARALLEL_LOOP(collapse=3) + do l = is3%beg, is3%end + do k = is2%beg, is2%end + do j = is1%beg, is1%end + $:GPU_LOOP(parallelism='[seq]') + do i = 1, sys_size flux_vf(i)%sf(j, k, l) = flux_${INFIX}$rsx_vf(j, k, l, i) end do end do diff --git a/src/simulation/m_weno.fpp b/src/simulation/m_weno.fpp index 89cc45f774..70e838379a 100644 --- a/src/simulation/m_weno.fpp +++ b/src/simulation/m_weno.fpp @@ -868,11 +868,12 @@ contains type(scalar_field), dimension(1:), intent(in) :: v_vf integer :: i, j, k, l, n_vars - $:GPU_PARALLEL_LOOP(collapse=4) - do i = 1, v_size - do l = idwbuff(3)%beg, idwbuff(3)%end - do k = idwbuff(2)%beg, idwbuff(2)%end - do j = idwbuff(1)%beg, idwbuff(1)%end + $:GPU_PARALLEL_LOOP(collapse=3) + do l = idwbuff(3)%beg, idwbuff(3)%end + do k = idwbuff(2)%beg, idwbuff(2)%end + do j = idwbuff(1)%beg, idwbuff(1)%end + $:GPU_LOOP(parallelism='[seq]') + do i = 1, v_size v_rs_weno(j, k, l, i) = v_vf(i)%sf(j, k, l) end do end do