From 33afdeca528becbe02d1f05604a730542b00f066 Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Mon, 12 Feb 2024 17:39:29 -0500 Subject: [PATCH] fix bug --- src/simulation/m_rhs.fpp | 187 ++++++++++++++++++++++----------------- 1 file changed, 108 insertions(+), 79 deletions(-) diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index c0142101a6..e8d0dff9df 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -872,9 +872,9 @@ contains do k = 0, m rhs_vf(j)%sf(k, l, q) = & rhs_vf(j)%sf(k, l, q) + 1d0/dx(k)* & - q_prim_vf(contxe + idir)%sf(k, l, q)* & - (flux_src_n_vf(j)%sf(k - 1, l, q) & - - flux_src_n_vf(j)%sf(k, l, q)) + q_prim_qp%vf(contxe + id)%sf(k, l, q)* & + (flux_src_n(1)%vf(j)%sf(k - 1, l, q) & + - flux_src_n(1)%vf(j)%sf(k, l, q)) end do end do end do @@ -889,9 +889,9 @@ contains do k = 0, m rhs_vf(j)%sf(k, l, q) = & rhs_vf(j)%sf(k, l, q) + 1d0/dx(k)* & - (q_cons_vf(j)%sf(k, l, q) - Kterm(k, l, q))* & - (flux_src_n_vf(j)%sf(k, l, q) & - - flux_src_n_vf(j)%sf(k - 1, l, q)) + (q_cons_qp%vf(j)%sf(k, l, q) - Kterm(k, l, q))* & + (flux_src_n(1)%vf(j)%sf(k, l, q) & + - flux_src_n(1)%vf(j)%sf(k - 1, l, q)) end do end do end do @@ -902,9 +902,9 @@ contains do k = 0, m rhs_vf(j)%sf(k, l, q) = & rhs_vf(j)%sf(k, l, q) + 1d0/dx(k)* & - (q_cons_vf(j)%sf(k, l, q) + Kterm(k, l, q))* & - (flux_src_n_vf(j)%sf(k, l, q) & - - flux_src_n_vf(j)%sf(k - 1, l, q)) + (q_cons_qp%vf(j)%sf(k, l, q) + Kterm(k, l, q))* & + (flux_src_n(1)%vf(j)%sf(k, l, q) & + - flux_src_n(1)%vf(j)%sf(k - 1, l, q)) end do end do end do @@ -918,9 +918,9 @@ contains do k = 0, m rhs_vf(j)%sf(k, l, q) = & rhs_vf(j)%sf(k, l, q) + 1d0/dx(k)* & - q_cons_vf(j)%sf(k, l, q)* & - (flux_src_n_vf(j)%sf(k, l, q) & - - flux_src_n_vf(j)%sf(k - 1, l, q)) + q_cons_qp%vf(j)%sf(k, l, q)* & + (flux_src_n(1)%vf(j)%sf(k, l, q) & + - flux_src_n(1)%vf(j)%sf(k - 1, l, q)) end do end do end do @@ -1016,15 +1016,15 @@ contains do q = 0, m rhs_vf(j)%sf(q, k, l) = & rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & - q_prim_vf(contxe + idir)%sf(q, k, l)* & - (flux_src_n_vf(j)%sf(q, k - 1, l) & - - flux_src_n_vf(j)%sf(q, k, l)) + q_prim_qp%vf(contxe + id)%sf(q, k, l)* & + (flux_src_n(2)%vf(j)%sf(q, k - 1, l) & + - flux_src_n(2)%vf(j)%sf(q, k, l)) end do end do end do end do else - + if (alt_soundspeed) then do j = advxb, advxe if ((j == advxe) .and. (bubbles .neqv. .true.)) then @@ -1034,9 +1034,9 @@ contains do q = 0, m rhs_vf(j)%sf(q, k, l) = & rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & - (q_cons_vf(j)%sf(q, k, l) - Kterm(q, k, l))* & - (flux_src_n_vf(j)%sf(q, k, l) & - - flux_src_n_vf(j)%sf(q, k - 1, l)) + (q_cons_qp%vf(j)%sf(q, k, l) - Kterm(q, k, l))* & + (flux_src_n(2)%vf(j)%sf(q, k, l) & + - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) end do end do end do @@ -1048,8 +1048,8 @@ contains rhs_vf(j)%sf(q, k, l) = & rhs_vf(j)%sf(q, k, l) - & (Kterm(q, k, l)/2d0/y_cc(k))* & - (flux_src_n_vf(j)%sf(q, k, l) & - + flux_src_n_vf(j)%sf(q, k - 1, l)) + (flux_src_n(2)%vf(j)%sf(q, k, l) & + + flux_src_n(2)%vf(j)%sf(q, k - 1, l)) end do end do end do @@ -1061,9 +1061,9 @@ contains do q = 0, m rhs_vf(j)%sf(q, k, l) = & rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & - (q_cons_vf(j)%sf(q, k, l) + Kterm(q, k, l))* & - (flux_src_n_vf(j)%sf(q, k, l) & - - flux_src_n_vf(j)%sf(q, k - 1, l)) + (q_cons_qp%vf(j)%sf(q, k, l) + Kterm(q, k, l))* & + (flux_src_n(2)%vf(j)%sf(q, k, l) & + - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) end do end do end do @@ -1075,8 +1075,8 @@ contains rhs_vf(j)%sf(q, k, l) = & rhs_vf(j)%sf(q, k, l) + & (Kterm(q, k, l)/2d0/y_cc(k))* & - (flux_src_n_vf(j)%sf(q, k, l) & - + flux_src_n_vf(j)%sf(q, k - 1, l)) + (flux_src_n(2)%vf(j)%sf(q, k, l) & + + flux_src_n(2)%vf(j)%sf(q, k - 1, l)) end do end do end do @@ -1091,9 +1091,9 @@ contains do q = 0, m rhs_vf(j)%sf(q, k, l) = & rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & - q_cons_vf(j)%sf(q, k, l)* & - (flux_src_n_vf(j)%sf(q, k, l) & - - flux_src_n_vf(j)%sf(q, k - 1, l)) + q_cons_qp%vf(j)%sf(q, k, l)* & + (flux_src_n(2)%vf(j)%sf(q, k, l) & + - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) end do end do end do @@ -1182,15 +1182,16 @@ contains if (grid_geometry == 3) then if (riemann_solver == 1) then + !$acc parallel loop collapse(4) gang vector default(present) do j = advxb, advxe - do k = 0, p - do q = 0, n - do l = 0, m - rhs_vf(j)%sf(l, q, k) = & - rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)/y_cc(q)* & - q_prim_vf(contxe + idir)%sf(l, q, k)* & - (flux_src_n_vf(j)%sf(l, q, k - 1) & - - flux_src_n_vf(j)%sf(l, q, k)) + do l = 0, p + do k = 0, n + do q = 0, m + rhs_vf(j)%sf(q, k, l) = & + rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & + q_prim_qp%vf(contxe + id)%sf(q, k, l)* & + (flux_src_n(2)%vf(j)%sf(q, k - 1, l) & + - flux_src_n(2)%vf(j)%sf(q, k, l)) end do end do end do @@ -1201,43 +1202,71 @@ contains do j = advxb, advxe if ((j == advxe) .and. (bubbles .neqv. .true.)) then !$acc parallel loop collapse(3) gang vector default(present) - do k = 0, p - do q = 0, n - do l = 0, m - rhs_vf(j)%sf(l, q, k) = & - rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)/y_cc(q)* & - (q_cons_vf(j)%sf(l, q, k) - Kterm(l, q, k))* & - (flux_src_n_vf(j)%sf(l, q, k) & - - flux_src_n_vf(j)%sf(l, q, k - 1)) + do l = 0, p + do k = 0, n + do q = 0, m + rhs_vf(j)%sf(q, k, l) = & + rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & + (q_cons_qp%vf(j)%sf(q, k, l) - Kterm(q, k, l))* & + (flux_src_n(2)%vf(j)%sf(q, k, l) & + - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) end do end do end do + if (cyl_coord) then + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do q = 0, m + rhs_vf(j)%sf(q, k, l) = & + rhs_vf(j)%sf(q, k, l) - & + (Kterm(q, k, l)/2d0/y_cc(k))* & + (flux_src_n(2)%vf(j)%sf(q, k, l) & + + flux_src_n(2)%vf(j)%sf(q, k - 1, l)) + end do + end do + end do + end if else if ((j == advxb) .and. (bubbles .neqv. .true.)) then !$acc parallel loop collapse(3) gang vector default(present) - do k = 0, p - do q = 0, n - do l = 0, m - rhs_vf(j)%sf(l, q, k) = & - rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)/y_cc(q)* & - (q_cons_vf(j)%sf(l, q, k) + Kterm(l, q, k))* & - (flux_src_n_vf(j)%sf(l, q, k) & - - flux_src_n_vf(j)%sf(l, q, k - 1)) + do l = 0, p + do k = 0, n + do q = 0, m + rhs_vf(j)%sf(q, k, l) = & + rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & + (q_cons_qp%vf(j)%sf(q, k, l) + Kterm(q, k, l))* & + (flux_src_n(2)%vf(j)%sf(q, k, l) & + - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) end do end do end do + if (cyl_coord) then + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do q = 0, m + rhs_vf(j)%sf(q, k, l) = & + rhs_vf(j)%sf(q, k, l) + & + (Kterm(q, k, l)/2d0/y_cc(k))* & + (flux_src_n(2)%vf(j)%sf(q, k, l) & + + flux_src_n(2)%vf(j)%sf(q, k - 1, l)) + end do + end do + end do + end if end if end do else !$acc parallel loop collapse(4) gang vector default(present) do j = advxb, advxe - do k = 0, p - do q = 0, n - do l = 0, m - rhs_vf(j)%sf(l, q, k) = & - rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)/y_cc(q)* & - q_cons_vf(j)%sf(l, q, k)* & - (flux_src_n_vf(j)%sf(l, q, k) & - - flux_src_n_vf(j)%sf(l, q, k - 1)) + do l = 0, p + do k = 0, n + do q = 0, m + rhs_vf(j)%sf(q, k, l) = & + rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & + q_cons_qp%vf(j)%sf(q, k, l)* & + (flux_src_n(2)%vf(j)%sf(q, k, l) & + - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) end do end do end do @@ -1253,15 +1282,15 @@ contains do l = 0, m rhs_vf(j)%sf(l, q, k) = & rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)* & - q_prim_vf(contxe + idir)%sf(l, q, k)* & - (flux_src_n_vf(j)%sf(l, q, k - 1) & - - flux_src_n_vf(j)%sf(l, q, k)) + q_prim_qp%vf(contxe + id)%sf(l, q, k)* & + (flux_src_n(3)%vf(j)%sf(l, q, k - 1) & + - flux_src_n(3)%vf(j)%sf(l, q, k)) end do end do end do end do else - + if (alt_soundspeed) then do j = advxb, advxe if ((j == advxe) .and. (bubbles .neqv. .true.)) then @@ -1271,9 +1300,9 @@ contains do l = 0, m rhs_vf(j)%sf(l, q, k) = & rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)* & - (q_cons_vf(j)%sf(l, q, k) - Kterm(l, q, k))* & - (flux_src_n_vf(j)%sf(l, q, k) & - - flux_src_n_vf(j)%sf(l, q, k - 1)) + (q_cons_qp%vf(j)%sf(l, q, k) - Kterm(l, q, k))* & + (flux_src_n(3)%vf(j)%sf(l, q, k) & + - flux_src_n(3)%vf(j)%sf(l, q, k - 1)) end do end do end do @@ -1284,9 +1313,9 @@ contains do l = 0, m rhs_vf(j)%sf(l, q, k) = & rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)* & - (q_cons_vf(j)%sf(l, q, k) + Kterm(l, q, k))* & - (flux_src_n_vf(j)%sf(l, q, k) & - - flux_src_n_vf(j)%sf(l, q, k - 1)) + (q_cons_qp%vf(j)%sf(l, q, k) + Kterm(l, q, k))* & + (flux_src_n(3)%vf(j)%sf(l, q, k) & + - flux_src_n(3)%vf(j)%sf(l, q, k - 1)) end do end do end do @@ -1300,9 +1329,9 @@ contains do l = 0, m rhs_vf(j)%sf(l, q, k) = & rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)* & - q_cons_vf(j)%sf(l, q, k)* & - (flux_src_n_vf(j)%sf(l, q, k) & - - flux_src_n_vf(j)%sf(l, q, k - 1)) + q_cons_qp%vf(j)%sf(l, q, k)* & + (flux_src_n(3)%vf(j)%sf(l, q, k) & + - flux_src_n(3)%vf(j)%sf(l, q, k - 1)) end do end do end do @@ -1429,21 +1458,21 @@ contains end subroutine s_compute_rhs ! ----------------------------------------- - ! subroutine s_compute_advection_source_term(idir, rhs_vf, q_cons_vf, q_prim_vf, flux_src_n_vf) + ! subroutine s_compute_advection_source_term(id, rhs_vf, q_cons_vf, q_prim_vf, flux_src_n_vf) ! type(scalar_field), dimension(sys_size), intent(INOUT) :: q_cons_vf ! type(scalar_field), dimension(sys_size), intent(INOUT) :: q_prim_vf ! type(scalar_field), dimension(sys_size), intent(INOUT) :: rhs_vf ! type(scalar_field), dimension(sys_size), intent(INOUT) :: flux_src_n_vf - ! integer :: idir + ! integer :: id ! integer :: i, j, k, l, q - ! if (idir == 1) then + ! if (id == 1) then - ! elseif (idir == 2) then + ! elseif (id == 2) then - ! elseif (idir == 3) then + ! elseif (id == 3) then ! end if