Skip to content

Commit

Permalink
5th order 3D TGV Works
Browse files Browse the repository at this point in the history
  • Loading branch information
@wilfonba
wilfonba committed Jan 20, 2025
1 parent f982447 commit 3bdcfba
Show file tree
Hide file tree
Showing 4 changed files with 106 additions and 104 deletions.
2 changes: 1 addition & 1 deletion examples/3D_TaylorGreenVortex/case.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@
"dt": dt,
"t_step_start": 0,
"t_step_stop": Nt,
"t_step_save": 1, #int(Nt / 100),
"t_step_save": int(Nt / 100),
# Simulation Algorithm Parameters
"num_patches": 1,
"model_eqns": 2,
Expand Down
2 changes: 1 addition & 1 deletion src/simulation/m_global_parameters.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1121,7 +1121,7 @@ contains
end if
if (igr) then
buff_size = max(buff_size, 4)
buff_size = 4
end if
if (probe_wrt) then
Expand Down
205 changes: 103 additions & 102 deletions src/simulation/m_igr.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ module m_igr
implicit none

private; public :: s_initialize_igr_module, &
s_reconstruct_deriv, &
s_reconstruct_igr, &
s_igr_jacobi_iteration, &
s_igr_riemann_solver, &
s_initialize_igr, &
Expand Down Expand Up @@ -100,7 +100,7 @@ contains

! WENO like reconstructions using only the optimial weights as given in
! Balsara & Shu (1999)
subroutine s_reconstruct_deriv(qL, qR, q_prim, idir)
subroutine s_reconstruct_igr(qL, qR, q_prim, idir)

real(wp), &
dimension(startx:, starty:, startz:), &
Expand Down Expand Up @@ -407,7 +407,7 @@ contains
end if
end if

end subroutine s_reconstruct_deriv
end subroutine s_reconstruct_igr

subroutine s_igr_jacobi_iteration()

Expand Down Expand Up @@ -961,6 +961,7 @@ contains
end do
end do


!$acc parallel loop collapse(3) gang vector default(present) private(rho_lx, rho_rx, rho_ly, rho_ry)
do l = 0, p
do k = idwbuff(2)%beg + 1, idwbuff(2)%end - 1
Expand All @@ -973,7 +974,6 @@ contains
fd_coeff(j,k,l) = 1._wp / rho_igr(j, k, l)

fd_coeff(j,k,l) = fd_coeff(j,k,l) + alf_igr * ( (1._wp / dx(j)**2._wp) * (rho_lx + rho_rx) + (1._wp / dy(k)**2._wp) *(rho_ly + rho_ry) )

end do
end do
end do
Expand Down Expand Up @@ -1001,7 +1001,6 @@ contains
8._wp*q_prim_vf(momxb)%sf(j,k,l-1) + &
q_prim_vf(momxb)%sf(j,k,l-2) - &
q_prim_vf(momxb)%sf(j,k,l+2) )

end do
end do
end do
Expand All @@ -1027,7 +1026,6 @@ contains
8._wp*q_prim_vf(momxb+1)%sf(j,k,l-1) + &
q_prim_vf(momxb+1)%sf(j,k,l-2) - &
q_prim_vf(momxb+1)%sf(j,k,l+2) )

end do
end do
end do
Expand All @@ -1053,7 +1051,6 @@ contains
8._wp*q_prim_vf(momxb+2)%sf(j,k,l-1) + &
q_prim_vf(momxb+2)%sf(j,k,l-2) - &
q_prim_vf(momxb+2)%sf(j,k,l+2) )

end do
end do
end do
Expand Down Expand Up @@ -1097,73 +1094,77 @@ contains
type(scalar_field), intent(in), dimension(sys_size) :: q_prim_vf
integer, intent(in) :: idir

if (idir == 1) then
if(p == 0) then
!$acc parallel loop collapse(4) gang vector default(present)
do i = 1, sys_size - 1
do l = 0, p
do k = idwbuff(2)%beg + 2, idwbuff(2)%end - 2
do j = idwbuff(1)%beg + 2, idwbuff(1)%end - 2
qL_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j-2, k, l) + 27._wp * q_prim_vf(i)%sf(j-1, k, l) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j+1, k, l) + 2._wp * q_prim_vf(i)%sf(j+2, k, l))
qR_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j+2, k, l) + 27._wp * q_prim_vf(i)%sf(j+1, k, l) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j-1, k, l) + 2._wp * q_prim_vf(i)%sf(j-2, k, l))
end do
end do
end do
end do
else
!$acc parallel loop collapse(4) gang vector default(present)
do i = 1, sys_size - 1
do l = idwbuff(3)%beg + 2, idwbuff(3)%end - 2
do k = idwbuff(2)%beg + 2, idwbuff(2)%end - 2
do j = idwbuff(1)%beg + 2, idwbuff(1)%end - 2
qL_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j-2, k, l) + 27._wp * q_prim_vf(i)%sf(j-1, k, l) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j+1, k, l) + 2._wp * q_prim_vf(i)%sf(j+2, k, l))
qR_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j+2, k, l) + 27._wp * q_prim_vf(i)%sf(j+1, k, l) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j-1, k, l) + 2._wp * q_prim_vf(i)%sf(j-2, k, l))
end do
end do
end do
end do
end if
else if (idir == 2) then
if(p == 0) then
!$acc parallel loop collapse(4) gang vector default(present)
do i = 1, sys_size - 1
do l = 0, p
do k = idwbuff(2)%beg + 2, idwbuff(2)%end - 2
do j = idwbuff(1)%beg + 2, idwbuff(1)%end - 2
qL_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j, k-2, l) + 27._wp * q_prim_vf(i)%sf(j, k-1, l) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j, k+1, l) + 2._wp * q_prim_vf(i)%sf(j, k+2, l))
qR_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j, k+2, l) + 27._wp * q_prim_vf(i)%sf(j, k+1, l) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j, k-1, l) + 2._wp * q_prim_vf(i)%sf(j, k-2, l))
end do
end do
end do
end do
else
!$acc parallel loop collapse(4) gang vector default(present)
do i = 1, sys_size - 1
do l = idwbuff(3)%beg + 2, idwbuff(3)%end - 2
do k = idwbuff(2)%beg + 2, idwbuff(2)%end - 2
do j = idwbuff(1)%beg + 2, idwbuff(1)%end - 2
qL_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j, k-2, l) + 27._wp * q_prim_vf(i)%sf(j, k-1, l) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j, k+1, l) + 2._wp * q_prim_vf(i)%sf(j, k+2, l))
qR_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j, k+2, l) + 27._wp * q_prim_vf(i)%sf(j, k+1, l) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j, k-1, l) + 2._wp * q_prim_vf(i)%sf(j, k-2, l))
end do
end do
end do
end do
end if
else if (idir == 3) then
!$acc parallel loop collapse(4) gang vector default(present)
do i = 1, sys_size - 1
do l = idwbuff(3)%beg + 2, idwbuff(3)%end - 2
do k = idwbuff(2)%beg + 2, idwbuff(2)%end - 2
do j = idwbuff(1)%beg + 2, idwbuff(1)%end - 2
qL_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j, k, l-2) + 27._wp * q_prim_vf(i)%sf(j, k, l-1) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j, k, l+1) + 2._wp * q_prim_vf(i)%sf(j, k, l+2))
qR_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j, k, l+2) + 27._wp * q_prim_vf(i)%sf(j, k, l+1) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j, k, l-1) + 2._wp * q_prim_vf(i)%sf(j, k, l-2))
end do
end do
end do
end do
end if
do i = 1, sys_size - 1
call s_reconstruct_igr(qL_rs_vf(:,:,:,i), qR_rs_vf(:,:,:,i), q_prim_vf(i)%sf, idir)
end do

call s_reconstruct_deriv(FL, FR, jac, idir)
!if (idir == 1) then
!if(p == 0) then
!!$acc parallel loop collapse(4) gang vector default(present)
!do i = 1, sys_size - 1
!do l = 0, p
!do k = idwbuff(2)%beg + 2, idwbuff(2)%end - 2
!do j = idwbuff(1)%beg + 2, idwbuff(1)%end - 2
!qL_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j-2, k, l) + 27._wp * q_prim_vf(i)%sf(j-1, k, l) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j+1, k, l) + 2._wp * q_prim_vf(i)%sf(j+2, k, l))
!qR_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j+2, k, l) + 27._wp * q_prim_vf(i)%sf(j+1, k, l) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j-1, k, l) + 2._wp * q_prim_vf(i)%sf(j-2, k, l))
!end do
!end do
!end do
!end do
!else
!!$acc parallel loop collapse(4) gang vector default(present)
!do i = 1, sys_size - 1
!do l = idwbuff(3)%beg + 2, idwbuff(3)%end - 2
!do k = idwbuff(2)%beg + 2, idwbuff(2)%end - 2
!do j = idwbuff(1)%beg + 2, idwbuff(1)%end - 2
!qL_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j-2, k, l) + 27._wp * q_prim_vf(i)%sf(j-1, k, l) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j+1, k, l) + 2._wp * q_prim_vf(i)%sf(j+2, k, l))
!qR_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j+2, k, l) + 27._wp * q_prim_vf(i)%sf(j+1, k, l) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j-1, k, l) + 2._wp * q_prim_vf(i)%sf(j-2, k, l))
!end do
!end do
!end do
!end do
!end if
!else if (idir == 2) then
!if(p == 0) then
!!$acc parallel loop collapse(4) gang vector default(present)
!do i = 1, sys_size - 1
!do l = 0, p
!do k = idwbuff(2)%beg + 2, idwbuff(2)%end - 2
!do j = idwbuff(1)%beg + 2, idwbuff(1)%end - 2
!qL_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j, k-2, l) + 27._wp * q_prim_vf(i)%sf(j, k-1, l) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j, k+1, l) + 2._wp * q_prim_vf(i)%sf(j, k+2, l))
!qR_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j, k+2, l) + 27._wp * q_prim_vf(i)%sf(j, k+1, l) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j, k-1, l) + 2._wp * q_prim_vf(i)%sf(j, k-2, l))
!end do
!end do
!end do
!end do
!else
!!$acc parallel loop collapse(4) gang vector default(present)
!do i = 1, sys_size - 1
!do l = idwbuff(3)%beg + 2, idwbuff(3)%end - 2
!do k = idwbuff(2)%beg + 2, idwbuff(2)%end - 2
!do j = idwbuff(1)%beg + 2, idwbuff(1)%end - 2
!qL_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j, k-2, l) + 27._wp * q_prim_vf(i)%sf(j, k-1, l) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j, k+1, l) + 2._wp * q_prim_vf(i)%sf(j, k+2, l))
!qR_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j, k+2, l) + 27._wp * q_prim_vf(i)%sf(j, k+1, l) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j, k-1, l) + 2._wp * q_prim_vf(i)%sf(j, k-2, l))
!end do
!end do
!end do
!end do
!end if
!else if (idir == 3) then
!!$acc parallel loop collapse(4) gang vector default(present)
!do i = 1, sys_size - 1
!do l = idwbuff(3)%beg + 2, idwbuff(3)%end - 2
!do k = idwbuff(2)%beg + 2, idwbuff(2)%end - 2
!do j = idwbuff(1)%beg + 2, idwbuff(1)%end - 2
!qL_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j, k, l-2) + 27._wp * q_prim_vf(i)%sf(j, k, l-1) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j, k, l+1) + 2._wp * q_prim_vf(i)%sf(j, k, l+2))
!qR_rs_vf(j, k, l, i) = (1._wp/60._wp) * (-3._wp * q_prim_vf(i)%sf(j, k, l+2) + 27._wp * q_prim_vf(i)%sf(j, k, l+1) + 47._wp * q_prim_vf(i)%sf(j, k, l) -13._wp * q_prim_vf(i)%sf(j, k, l-1) + 2._wp * q_prim_vf(i)%sf(j, k, l-2))
!end do
!end do
!end do
!end do
!end if

call s_reconstruct_igr(FL, FR, jac, idir)

end subroutine s_reconstruct_prim_vars_igr

Expand All @@ -1172,41 +1173,41 @@ contains
integer, intent(in) :: idir

if (idir == 1) then
call s_reconstruct_deriv(duLx, duRx, dux, 1)
call s_reconstruct_deriv(dvLx, dvRx, dvx, 1)
call s_reconstruct_deriv(dwLx, dwRx, dwx, 1)
call s_reconstruct_igr(duLx, duRx, dux, 1)
call s_reconstruct_igr(dvLx, dvRx, dvx, 1)
call s_reconstruct_igr(dwLx, dwRx, dwx, 1)

call s_reconstruct_deriv(duLy, duRy, duy, 1)
call s_reconstruct_deriv(dvLy, dvRy, dvy, 1)
call s_reconstruct_deriv(dwLy, dwRy, dwy, 1)
call s_reconstruct_igr(duLy, duRy, duy, 1)
call s_reconstruct_igr(dvLy, dvRy, dvy, 1)
call s_reconstruct_igr(dwLy, dwRy, dwy, 1)

call s_reconstruct_deriv(duLz, duRz, duz, 1)
call s_reconstruct_deriv(dvLz, dvRz, dvz, 1)
call s_reconstruct_deriv(dwLz, dwRz, dwz, 1)
call s_reconstruct_igr(duLz, duRz, duz, 1)
call s_reconstruct_igr(dvLz, dvRz, dvz, 1)
call s_reconstruct_igr(dwLz, dwRz, dwz, 1)
elseif (idir == 2) then
call s_reconstruct_deriv(duLx, duRx, dux, 2)
call s_reconstruct_deriv(dvLx, dvRx, dvx, 2)
call s_reconstruct_deriv(dwLx, dwRx, dwx, 2)
call s_reconstruct_igr(duLx, duRx, dux, 2)
call s_reconstruct_igr(dvLx, dvRx, dvx, 2)
call s_reconstruct_igr(dwLx, dwRx, dwx, 2)

call s_reconstruct_deriv(duLy, duRy, duy, 2)
call s_reconstruct_deriv(dvLy, dvRy, dvy, 2)
call s_reconstruct_deriv(dwLy, dwRy, dwy, 2)
call s_reconstruct_igr(duLy, duRy, duy, 2)
call s_reconstruct_igr(dvLy, dvRy, dvy, 2)
call s_reconstruct_igr(dwLy, dwRy, dwy, 2)

call s_reconstruct_deriv(duLz, duRz, duz, 2)
call s_reconstruct_deriv(dvLz, dvRz, dvz, 2)
call s_reconstruct_deriv(dwLz, dwRz, dwz, 2)
call s_reconstruct_igr(duLz, duRz, duz, 2)
call s_reconstruct_igr(dvLz, dvRz, dvz, 2)
call s_reconstruct_igr(dwLz, dwRz, dwz, 2)
elseif (idir == 3) then
call s_reconstruct_deriv(duLx, duRx, dux, 3)
call s_reconstruct_deriv(dvLx, dvRx, dvx, 3)
call s_reconstruct_deriv(dwLx, dwRx, dwx, 3)
call s_reconstruct_igr(duLx, duRx, dux, 3)
call s_reconstruct_igr(dvLx, dvRx, dvx, 3)
call s_reconstruct_igr(dwLx, dwRx, dwx, 3)

call s_reconstruct_deriv(duLy, duRy, duy, 3)
call s_reconstruct_deriv(dvLy, dvRy, dvy, 3)
call s_reconstruct_deriv(dwLy, dwRy, dwy, 3)
call s_reconstruct_igr(duLy, duRy, duy, 3)
call s_reconstruct_igr(dvLy, dvRy, dvy, 3)
call s_reconstruct_igr(dwLy, dwRy, dwy, 3)

call s_reconstruct_deriv(duLz, duRz, duz, 3)
call s_reconstruct_deriv(dvLz, dvRz, dvz, 3)
call s_reconstruct_deriv(dwLz, dwRz, dwz, 3)
call s_reconstruct_igr(duLz, duRz, duz, 3)
call s_reconstruct_igr(dvLz, dvRz, dvz, 3)
call s_reconstruct_igr(dwLz, dwRz, dwz, 3)
end if

end subroutine s_get_viscous_igr
Expand Down Expand Up @@ -1240,8 +1241,8 @@ contains
do j = 0, m
rhs_vf(i)%sf(j, k, l) = &
rhs_vf(i)%sf(j, k, l) + 1._wp/dy(k)* &
(flux_vf(i)%sf(k - 1, j , l) &
- flux_vf(i)%sf(k, j , l))
(flux_vf(i)%sf(j, k - 1, l) &
- flux_vf(i)%sf(j, k , l))
end do
end do
end do
Expand Down
1 change: 1 addition & 0 deletions toolchain/mfc/case.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -212,6 +212,7 @@ def __get_sim_fpp(self, print: bool) -> str:
#:set wenoz = {wenoz}
#:set teno = {teno}
#:set wenoz_q = {self.params.get("wenoz_q", -1)}
#:set igr_order = {int(self.params.get("igr_order", 1))}
"""

return """\
Expand Down

0 comments on commit 3bdcfba

Please sign in to comment.