From 23d00a9e8e4fd7c855a48918de9fc4bab24141c9 Mon Sep 17 00:00:00 2001 From: Anand Date: Wed, 17 Apr 2024 17:43:37 -0400 Subject: [PATCH] NEW LW --- src/common/m_mpi_common.fpp | 2 + src/post_process/m_start_up.f90 | 29 +- src/pre_process/m_data_output.fpp | 4 +- src/simulation/m_rhs.fpp | 480 ++++++----------------------- src/simulation/m_time_steppers.fpp | 6 +- 5 files changed, 113 insertions(+), 408 deletions(-) diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index 4d3064d90..77ed2d48e 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -110,6 +110,7 @@ contains end if end if + ! Define the view for each variable do i = 1, sys_size call MPI_TYPE_CREATE_SUBARRAY(num_dims, sizes_glb, sizes_loc, start_idx, & @@ -128,6 +129,7 @@ contains end if #endif + #endif end subroutine s_initialize_mpi_data ! --------------------------------- diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90 index a136d8387..8644905df 100644 --- a/src/post_process/m_start_up.f90 +++ b/src/post_process/m_start_up.f90 @@ -289,20 +289,21 @@ subroutine s_save_data(t_step, varname, pres, c, H) -offset_z%beg:p + offset_z%end) q_sf = 2.5d0*q_prim_vf(E_idx)%sf + 0.5d0*q_prim_vf(1)%sf(j,k,l)*& - (q_prim_vf(mom_idx%beg)%sf**2d0 + q_prim_vf(mom_idx%beg + 1)%sf**2d0) - - En_tot = 0d0 - rho_tot = 0d0 - do k = 0, n - do j = 0, m - En_tot = En_tot + q_sf(j,k,0) - rho_tot = rho_tot + q_cons_vf(1)%sf(j,k,0) - end do - end do - En_tot = En_tot/(m+1)**2d0 - rho_tot = 0.5d0*(rho_tot/(m+1)**2d0)**1.4 - print *, "POT AVG", rho_tot - print *, "En_tot", En_tot + (q_prim_vf(mom_idx%beg)%sf**2d0 + q_prim_vf(mom_idx%beg + 1)%sf**2d0) + + En_tot = 0d0 + rho_tot = 0d0 + do k = 0, n + do j = 0, m + En_tot = En_tot + q_sf(j,k,0) + rho_tot = rho_tot + q_cons_vf(1)%sf(j,k,0) + end do + end do + En_tot = En_tot/((m+1)*(n+1)) + !En_tot = En_tot - 0.59012242316857255 + rho_tot = 0.5d0*(rho_tot/((m+1)*(n+1)))**1.4 + print *, "POT AVG", rho_tot + print *, "En_tot", En_tot inquire (FILE='En_tot.dat', EXIST=file_exists) if (file_exists) then open (1, file='En_tot.dat', position='append', status='old') diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp index 9c57d9d76..0493a681d 100644 --- a/src/pre_process/m_data_output.fpp +++ b/src/pre_process/m_data_output.fpp @@ -412,7 +412,9 @@ contains call MPI_FILE_DELETE(file_loc, mpi_info_int, ierr) end if call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), & - mpi_info_int, ifile, ierr) + mpi_info_int, ifile, ierr) + + ! Size of local arrays data_size = (m + 1)*(n + 1)*(p + 1) diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index 58a7d364f..14a5380ca 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -581,9 +581,9 @@ contains do l = iy%beg, iy%end do j = ix%beg, ix%end do i = 1, 2 - jac_igr(i, j, k) = 0d0 - jac_old_igr(i, j, k) = 0d0 - jac_rhs_igr(i, j, k) = 0d0 + jac_igr(i, j, 0) = 0d0 + jac_old_igr(i, j, 0) = 0d0 + jac_rhs_igr(i, j, 0) = 0d0 end do end do end do @@ -779,11 +779,7 @@ contains bcxb = bc_x%beg; bcxe = bc_x%end; bcyb = bc_y%beg; bcye = bc_y%end -! bcxb = -1; bcxe = -1; bcyb = -1; bcye = -1 - - !print *, 'BC', bcxb, bcxe, bcyb, bcye - - !bcxb = -1; bcxe = -1 + bcxb = -1; bcxe = -1; bcyb = -1; bcye = -1 !$acc update device(ix, iy, iz, bcxb, bcxe, bcyb, bcye) @@ -971,21 +967,13 @@ contains end do end do - - ! !$acc update host(dux_igr, duy_igr, dvx_igr, dvy_igr, F_igr) - ! print *, maxval(abs(dux_igr(0:m,0:n))), maxval(abs(duy_igr(0:m,0:n))), maxval(abs(dvx_igr(0:m,0:n))),maxval(abs(dvy_igr(0:m,0:n))) - ! print *, maxval(abs(F_igr(1, -1:m+1, -1:n+1))), maxval(abs(F_igr(2, -1:m+1, -1:n+1))), maxval(abs(F_igr(3, -1:m+1, -1:n+1))), maxval(abs(F_igr(4, -1:m+1, -1:n+1))) - - ! !$acc update host(rho_igr, alf_igr) - ! print *, maxval(abs(rho_igr)) - ! print *, alf_igr !$acc parallel loop collapse(2) gang vector default(present) do k = iy%beg + 1, iy%end - 1 do j = ix%beg + 1, ix%end - 1 - rhs_igr(1, j, k) = rho_igr(j, k) *alf_igr *2d0 * (dux_igr(j, k)**2d0 + duy_igr(j, k) * dvx_igr(j, k)) - rhs_igr(2, j, k) = rho_igr(j, k) *alf_igr *2d0 * (dux_igr(j, k)*duy_igr(j, k) + duy_igr(j, k) * dvy_igr(j, k)) - rhs_igr(3, j, k) = rho_igr(j, k) *alf_igr *2d0 * (dvx_igr(j, k)*dux_igr(j, k) + dvy_igr(j, k) * dvx_igr(j, k)) - rhs_igr(4, j, k) = rho_igr(j, k) *alf_igr *2d0 * (dvx_igr(j, k)*duy_igr(j, k) + dvy_igr(j, k)**2d0) + rhs_igr(1, j, k) = (dux_igr(j, k)**2d0 + duy_igr(j, k) * dvx_igr(j, k)) + rhs_igr(2, j, k) = (dux_igr(j, k)*duy_igr(j, k) + duy_igr(j, k) * dvy_igr(j, k)) + rhs_igr(3, j, k) = (dvx_igr(j, k)*dux_igr(j, k) + dvy_igr(j, k) * dvx_igr(j, k)) + rhs_igr(4, j, k) = (dvx_igr(j, k)*duy_igr(j, k) + dvy_igr(j, k)**2d0) end do end do @@ -993,30 +981,26 @@ contains ! print *, "MAX", maxval(abs(rhs_igr(1:4, 0:m,0:n))) !$acc parallel loop collapse(2) gang vector default(present) - do k = iy%beg + 2, iy%end - 1 - do j = ix%beg + 2, ix%end - 1 - rho_lx = 1d0 * (rho_igr(j,k) + rho_igr(j-1,k)) / 2d0 - rho_ly = 1d0 * (rho_igr(j,k) + rho_igr(j,k-1)) / 2d0 - - jac_rhs_igr(1, j, k) = (1d0 / dx(j))*(rhs_igr(1,j-1, k) - rhs_igr(1,j,k)) + (1d0 / dy(k))*(rhs_igr(2,j, k-1) - rhs_igr(2,j,k)) - jac_rhs_igr(2, j, k) = (1d0 / dx(j))*(rhs_igr(3,j-1, k) - rhs_igr(3,j,k)) + (1d0 / dy(k))*(rhs_igr(4,j, k-1) - rhs_igr(4,j,k)) + do k = iy%beg + 1, iy%end - 1 + do j = ix%beg + 1, ix%end - 1 + jac_rhs_igr(1, j, k) = alf_igr * (rhs_igr(1, j, k) + rhs_igr(4, j, k) + (dux_igr(j, k) + dvy_igr(j, k))**2d0 ) end do end do !$acc parallel loop collapse(2) gang vector default(present) private(rho_lx, rho_rx, rho_ly, rho_ry) do k = 0, n do j = 0, m - rho_lx = 1d0 * (rho_igr(j,k) + rho_igr(j-1,k)) / 2d0 - rho_rx = 1d0 * (rho_igr(j,k) + rho_igr(j+1,k)) / 2d0 - rho_ly = 1d0 * (rho_igr(j,k) + rho_igr(j,k-1)) / 2d0 - rho_ry = 1d0 * (rho_igr(j,k) + rho_igr(j,k+1)) / 2d0 - - rho_lx = rho_igr(j-1,k) - rho_rx = rho_igr(j,k) - rho_ly = rho_igr(j,k-1) - rho_ry = rho_igr(j,k) - - fd_coeff(j,k) = rho_igr(j,k) / alf_igr + rho_lx = 0.5d0 *(1d0 / rho_igr(j,k) + 1d0 / rho_igr(j-1,k)) + rho_rx = 0.5d0 *(1d0 / rho_igr(j,k) + 1d0 / rho_igr(j+1,k)) + rho_ly = 0.5d0 *(1d0 / rho_igr(j,k) + 1d0 / rho_igr(j,k-1)) + rho_ry = 0.5d0 *(1d0 / rho_igr(j,k) + 1d0 / rho_igr(j,k+1)) + + ! rho_lx = 1d0 / rho_igr(j-1, k) + ! rho_rx = 1d0 / rho_igr(j, k) + ! rho_ly = 1d0 / rho_igr(j, k-1) + ! rho_ry = 1d0 / rho_igr(j, k) + + fd_coeff(j,k) = 1d0 / rho_igr(j, k) if(bcxb < -1 .and. j == 0) then rho_lx = 0d0 @@ -1035,33 +1019,30 @@ contains rho_ry = 0d0 end if - fd_coeff(j,k) = fd_coeff(j,k) + (1d0 / dx(j)**2d0) * (rho_lx + rho_rx) + (1d0 / dy(k)**2d0) *(rho_ly + rho_ry) + fd_coeff(j,k) = fd_coeff(j,k) + alf_igr * ( (1d0 / dx(j)**2d0) * (rho_lx + rho_rx) + (1d0 / dy(k)**2d0) *(rho_ly + rho_ry) ) end do end do - !print *, "PR INIT", proc_rank, jac_old_igr(1, 0:buff_size-1, n), jac_old_igr(1, m+1:m+buff_size, n) - - omega = 0.67 + omega = 1d0 !$acc update device(omega) - - !print *, "RHO DIFF", maxval(abs(rho_igr(100, 0:n+1) - rho_igr(100, -1:n))) do i = 1, 10 !$acc parallel loop gang vector collapse(3) default(present) private(rho_lx, rho_ly, rho_rx, rho_ry) - do q = 1, 2 + do q = 1, 1 do k = 0, n do j = 0, m - rho_lx = 1d0 * (rho_igr(j,k) + rho_igr(j-1,k)) / 2d0 - rho_rx = 1d0 * (rho_igr(j,k) + rho_igr(j+1,k)) / 2d0 - rho_ly = 1d0 * (rho_igr(j,k) + rho_igr(j,k-1)) / 2d0 - rho_ry = 1d0 * (rho_igr(j,k) + rho_igr(j,k+1)) / 2d0 + rho_lx = 0.5d0 *(1d0 / rho_igr(j,k) + 1d0 / rho_igr(j-1,k)) + rho_rx = 0.5d0 *(1d0 / rho_igr(j,k) + 1d0 / rho_igr(j+1,k)) + rho_ly = 0.5d0 *(1d0 / rho_igr(j,k) + 1d0 / rho_igr(j,k-1)) + rho_ry = 0.5d0 *(1d0 / rho_igr(j,k) + 1d0 / rho_igr(j,k+1)) + + ! rho_lx = 1d0 / rho_igr(j-1, k) + ! rho_rx = 1d0 / rho_igr(j, k) + ! rho_ly = 1d0 / rho_igr(j, k-1) + ! rho_ry = 1d0 / rho_igr(j, k) - rho_lx = rho_igr(j-1,k) - rho_rx = rho_igr(j,k) - rho_ly = rho_igr(j,k-1) - rho_ry = rho_igr(j,k) jac_igr(q, j, k) = jac_rhs_igr(q, j, k) @@ -1082,8 +1063,8 @@ contains rho_ry = 0d0 end if - jac_igr(q,j,k) = jac_igr(q,j,k)+ (1d0 / dx(j)**2d0) * (rho_lx* jac_old_igr(q,j-1,k) + rho_rx*jac_old_igr(q,j+1,k)) & - + (1d0 / dy(k)**2d0) * (rho_ly* jac_old_igr(q,j,k-1) + rho_ry* jac_old_igr(q,j,k+1)) + jac_igr(q,j,k) = jac_igr(q,j,k) + alf_igr * (1d0 / dx(j)**2d0) * (rho_lx* jac_old_igr(q,j-1,k) + rho_rx*jac_old_igr(q,j+1,k)) & + + alf_igr * (1d0 / dy(k)**2d0) * (rho_ly* jac_old_igr(q,j,k-1) + rho_ry* jac_old_igr(q,j,k+1)) jac_igr(q, j, k) = omega * (1 / fd_coeff(j,k))*jac_igr(q,j,k) + (1 - omega)*jac_old_igr(q, j, k) @@ -1092,17 +1073,12 @@ contains end do end do - !call MPI_Barrier(MPI_COMM_WORLD) - - !print *, "PR INIT", proc_rank, fd_coeff(0:buff_size-1, 199), fd_coeff(m+1:m+buff_size, 199) - - if(bcxb >= -1) then if(bcxb >= 0) then - !call s_mpi_sendrecv_F_igr(jac_igr, 1, -1) + call s_mpi_sendrecv_F_igr(jac_igr, 1, -1) else !$acc parallel loop gang vector collapse(3) default(present) - do q = 1, 2 + do q = 1, 1 do k = 0, n do j = 1, buff_size jac_igr(q,-j, k) = jac_igr(q,m-j+1,k) @@ -1115,10 +1091,10 @@ contains if(bcxe >= -1) then if(bcxe >= 0) then - !call s_mpi_sendrecv_F_igr(jac_igr, 1, 1) + call s_mpi_sendrecv_F_igr(jac_igr, 1, 1) else !$acc parallel loop gang vector collapse(3) default(present) - do q = 1, 2 + do q = 1, 1 do k = 0, n do j = 1, buff_size jac_igr(q,m+j, k) = jac_igr(q,j-1,k) @@ -1130,10 +1106,10 @@ contains if(bcyb >= -1) then if(bcyb >= 0) then - !call s_mpi_sendrecv_F_igr(jac_igr, 2, -1) + call s_mpi_sendrecv_F_igr(jac_igr, 2, -1) else !$acc parallel loop gang vector collapse(3) default(present) - do q = 1, 2 + do q = 1, 1 do k = 1, buff_size do j = ix%beg, ix%end jac_igr(q,j,-k) = jac_igr(q,j,n-k+1) @@ -1145,10 +1121,10 @@ contains if(bcye >= -1) then if(bcye >= 0) then - !call s_mpi_sendrecv_F_igr(jac_igr, 2, 1) + call s_mpi_sendrecv_F_igr(jac_igr, 2, 1) else !$acc parallel loop gang vector collapse(3) default(present) - do q = 1, 2 + do q = 1, 1 do k = 1, buff_size do j = ix%beg, ix%end jac_igr(q,j,n+k) = jac_igr(q,j,k-1) @@ -1158,37 +1134,20 @@ contains end if end if - !if(i == 25) then - !!$acc update host(jac_igr, jac_old_igr, fd_coeff) - ! q = 1 - ! print *, "ITER",q, maxval(abs(jac_igr(q, -1:m+1, -1:n+1) - jac_old_igr(q,-1:m+1, -1:n+1))), maxval(abs( jac_old_igr(q,-1:m+1,-1:n+1))), maxval(abs(fd_coeff(-1:m+1, -1:n+1))) - ! q = 2 - ! print *, "ITER",q, maxval(abs(jac_igr(q, -1:m+1, -1:n+1) - jac_old_igr(q,-1:m+1, -1:n+1))), maxval(abs( jac_old_igr(q,-1:m+1,-1:n+1))), maxval(abs(fd_coeff(-1:m+1, -1:n+1))) - - !end if - !print *, "PR INIT", proc_rank, jac_old_igr(q, 0:buff_size-1, 199), jac_old_igr(q, m+1:m+buff_size, 199) - - !$acc parallel loop gang vector collapse(3) default(present) - do q = 1, 2 + do q = 1, 1 do k = iy%beg, iy%end do j = ix%beg, ix%end jac_old_igr(q, j, k) = jac_igr(q, j, k) end do end do - end do - - !print *, "PR", proc_rank, jac_old_igr(q, -buff_size:-1, 98), jac_igr(q, m-buff_size+1:m, 98) - !print *, "PR", proc_rank, jac_old_igr(q, 0:buff_size-1, 199), jac_old_igr(q, m+1:m+buff_size, 199) - !print *, "NAN", proc_rank, fd_coeff(0:buff_size-1, 199), fd_coeff(m:m+buff_size, 199) - !print *, "NAN", proc_rank, rhs_igr(3, m-buff_size+1:m, 98), rhs_igr(4, m-buff_size+1:m, 98) - + end do end do if(bcxb < -1) then !$acc parallel loop gang vector collapse(3) default(present) - do q = 1, 2 + do q = 1, 1 do k = 0, n do j = 1, buff_size jac_igr(q,-j,k) = jac_igr(q,0,k) @@ -1199,7 +1158,7 @@ contains if(bcxe < -1) then !$acc parallel loop gang vector collapse(3) default(present) - do q = 1, 2 + do q = 1, 1 do k = 0, n do j = 1, buff_size jac_igr(q,m+j,k) = jac_igr(q,m,k) @@ -1210,7 +1169,7 @@ contains if(bcyb < -1) then !$acc parallel loop gang vector collapse(3) default(present) - do q = 1, 2 + do q = 1, 1 do k = 1, buff_size do j = ix%beg, ix%end jac_igr(q,j,-k) = jac_igr(q,j,0) @@ -1221,7 +1180,7 @@ contains if(bcye < -1) then !$acc parallel loop gang vector collapse(3) default(present) - do q = 1, 2 + do q = 1, 1 do k = 1, buff_size do j = ix%beg, ix%end jac_igr(q,j,n+k) = jac_igr(q,j,n) @@ -1231,20 +1190,14 @@ contains end if !$acc parallel loop gang vector collapse(2) default(present) - do k = iy%beg , iy%end - 1 - do j = ix%beg , ix%end - 1 - F_igr(1, j, k) = 1d0* (rhs_igr(1, j, k) - rho_igr(j,k)*(1d0 / dx(j))*(jac_igr(1,j+1,k) - jac_igr(1,j,k))) - F_igr(2, j, k) = 1d0* (rhs_igr(2, j, k) - rho_igr(j,k)*(1d0 / dy(k))*(jac_igr(1,j,k+1) - jac_igr(1,j,k))) - F_igr(3, j, k) = 1d0* (rhs_igr(3, j, k) - rho_igr(j,k)*(1d0 / dx(j))*(jac_igr(2,j+1,k) - jac_igr(2,j,k))) - F_igr(4, j, k) = 1d0* (rhs_igr(4, j, k) - rho_igr(j,k)*(1d0 / dy(k))*(jac_igr(2,j,k+1) - jac_igr(2,j,k))) + do k = iy%beg , iy%end + do j = ix%beg , ix%end + F_igr(1, j, k) = jac_igr(1, j, k) end do end do - ! !$acc update host(F_igr) - ! print *, "F", maxval(abs(F_igr(1:4, -1:m+1, -1:n+1))) - !$acc parallel loop gang vector collapse(2) default(present) - do k = 0, n + do k = -1, n+1 do j = -1, m+1 flux_n(id)%vf(1)%sf(j,k,0) = & @@ -1259,26 +1212,30 @@ contains F_igr(1, j, k) flux_n(id)%vf(3)%sf(j, k, 0) = q_prim_qp%vf(contxb)%sf(j,k,0) * & - q_prim_qp%vf( momxb)%sf(j,k,0)*q_prim_qp%vf( momxb + 1)%sf(j,k,0) + & - F_igr(3, j, k) + q_prim_qp%vf( momxb)%sf(j,k,0)*q_prim_qp%vf( momxb + 1)%sf(j,k,0) end do end do - ! do i = 1, 3 - ! !$acc update host(flux_n(1)%vf(i)%sf) - ! end do - ! print *, "Flux rho X", maxval(abs(flux_n(1)%vf(1)%sf( -1:m+1, -1:n+1, 0))) - ! print *, "Flux u X", maxval(abs(flux_n(1)%vf(2)%sf( -1:m+1, -1:n+1, 0))) - ! print *, "Flux v X", maxval(abs(flux_n(1)%vf(3)%sf( -1:m+1, -1:n+1, 0))) - !$acc parallel loop collapse(3) gang vector default(present) do i = 1, 3 do k = 0, n do j = 0, m - rhs_vf(i)%sf(j,k,0) = & - 1d0/(2d0*dx(j)) * & - ( flux_n(id)%vf(i)%sf(j-1,k,0) - & - flux_n(id)%vf(i)%sf(j+1,k,0) ) + + if(lw == 1) then + rhs_vf(i)%sf(j,k,0) = & + 1d0/(2d0*dx(j)) * & + ( flux_n(id)%vf(i)%sf(j,k,0) - & + flux_n(id)%vf(i)%sf(j+1,k,0) + & + flux_n(id)%vf(i)%sf(j,k+1, 0) - & + flux_n(id)%vf(i)%sf(j+1,k+1, 0)) + else if(lw == 2) then + rhs_vf(i)%sf(j,k,0) = & + 1d0/(2d0*dx(j)) * & + ( flux_n(id)%vf(i)%sf(j-1,k-1,0) - & + flux_n(id)%vf(i)%sf(j,k-1,0) + & + flux_n(id)%vf(i)%sf(j-1,k, 0) - & + flux_n(id)%vf(i)%sf(j,k, 0)) + end if if(bcxb < -1) then rhs_vf(i)%sf(0,k,0) = 0d0 @@ -1299,13 +1256,6 @@ contains end do end do end do - ! do i = 1, 3 - ! !$acc update host(rhs_vf(i)%sf) - ! end do - - ! print *, "RHS RHO X", maxval(abs(rhs_vf(1)%sf(0:m, 0:n, 0))) - ! print *, "RHS U X", maxval(abs(rhs_vf(2)%sf(0:m, 0:n, 0))) - ! print *, "RHS V X", maxval(abs(rhs_vf(3)%sf(0:m, 0:n, 0))) end if if(.not. barotropic) then @@ -1323,215 +1273,6 @@ contains end do end if - !!Non-polytropic qbmm needs to account for change in bubble radius due to a change in nb - !if(qbmm .and. (.not. polytropic) ) then - ! !$acc parallel loop collapse(5) gang vector default(present) private(nb_q, nR, nR2, R, R2, nb_dot, nR_dot, nR2_dot, var) - ! do i = 1, nb - ! do q = 1, nnode - ! do l = 0, p - ! do k = 0, n - ! do j = 0, m - ! nb_q = q_cons_qp%vf(bubxb + (i-1)*nmom)%sf(j, k, l) - ! nR = q_cons_qp%vf(bubxb + 1 + (i-1)*nmom)%sf(j, k, l) - ! nR2 = q_cons_qp%vf(bubxb + 3 + (i-1)*nmom)%sf(j, k, l) - - ! R = q_prim_qp%vf(bubxb + 1 + (i-1)*nmom)%sf(j, k, l) - ! R2 = q_prim_qp%vf(bubxb + 3 + (i-1)*nmom)%sf(j, k, l) - - ! nb_dot = flux_n(1)%vf(bubxb + (i-1)*nmom)%sf(j - 1, k, l) - flux_n(1)%vf(bubxb + (i-1)*nmom)%sf(j , k, l) - ! nR_dot = flux_n(1)%vf(bubxb + 1 + (i-1)*nmom)%sf(j - 1, k, l) - flux_n(1)%vf(bubxb + 1 + (i-1)*nmom)%sf(j , k, l) - ! nR2_dot = flux_n(1)%vf(bubxb + 3 + (i-1)*nmom)%sf(j - 1, k, l) - flux_n(1)%vf(bubxb + 3 + (i-1)*nmom)%sf(j , k, l) - - ! rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0 * gam/ (dx(j) * R * nb_q ** 2 )* & - ! (nR_dot * nb_q - nR * nb_dot) * (pb(j, k, l, q, i)) - - ! if(R2 - R**2d0 > 0d0) then - ! var = R2 - R**2d0 - ! else - ! var = verysmall - ! end if - - ! if(q <= 2) then - ! rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0 * gam/ (dx(j) * R * nb_q ** 2 * dsqrt(var)) * & - ! (nR2_dot * nb_q - nR2 * nb_dot ) * (pb(j, k, l, q, i)) - ! rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) + 3d0 * gam/ (dx(j) * R * nb_q ** 2 * dsqrt(var)) * & - ! ( - 2d0 * (nR / nb_q) * (nR_dot * nb_q - nR * nb_dot )) * (pb(j, k, l, q, i)) - - ! else - ! rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0 * gam/ (dx(j) * R * nb_q ** 2 * dsqrt(var)) * & - ! (nR2_dot * nb_q - nR2 * nb_dot ) * (pb(j, k, l, q, i)) - ! rhs_pb(j, k, l, q, i) = rhs_pb(j, k, l, q, i) - 3d0 * gam/ (dx(j) * R * nb_q ** 2 * dsqrt(var)) * & - ! ( - 2d0 * (nR / nb_q) * (nR_dot * nb_q - nR * nb_dot )) * (pb(j, k, l, q, i)) - ! end if - - ! end do - ! end do - ! end do - ! end do - ! end do - !end if - - !if (riemann_solver == 1) then - ! !$acc parallel loop collapse(4) gang vector default(present) - ! do j = advxb, advxe - ! do q = 0, p - ! do l = 0, n - ! do k = 0, m - ! rhs_vf(j)%sf(k, l, q) = & - ! rhs_vf(j)%sf(k, l, q) + 1d0/dx(k)* & - ! 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 - ! end do - !else - ! if (alt_soundspeed) then - ! do j = advxb, advxe - ! if ((j == advxe) .and. (bubbles .neqv. .true.)) then - ! !$acc parallel loop collapse(3) gang vector default(present) - ! do q = 0, p - ! do l = 0, n - ! do k = 0, m - ! rhs_vf(j)%sf(k, l, q) = & - ! rhs_vf(j)%sf(k, l, q) + 1d0/dx(k)* & - ! (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 - ! else if ((j == advxb) .and. (bubbles .neqv. .true.)) then - ! !$acc parallel loop collapse(3) gang vector default(present) - ! do q = 0, p - ! do l = 0, n - ! do k = 0, m - ! rhs_vf(j)%sf(k, l, q) = & - ! rhs_vf(j)%sf(k, l, q) + 1d0/dx(k)* & - ! (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 - ! end if - ! end do - ! else - ! !$acc parallel loop collapse(4) gang vector default(present) - ! do j = advxb, advxe - ! do q = 0, p - ! do l = 0, n - ! do k = 0, m - ! rhs_vf(j)%sf(k, l, q) = & - ! rhs_vf(j)%sf(k, l, q) + 1d0/dx(k)* & - ! 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 - ! end do - ! end if - !end if - - !if (bubbles) then - ! if (qbmm) then - - ! !$acc parallel loop collapse(3) gang vector default(present) - ! do l = 0, p - ! do q = 0, n - ! do i = 0, m - - ! rhs_vf(alf_idx)%sf(i, q, l) = rhs_vf(alf_idx)%sf(i, q, l) + mom_sp(2)%sf(i, q, l) - ! j = bubxb - ! !$acc loop seq - ! do k = 1, nb - ! rhs_vf(j)%sf(i, q, l) = & - ! rhs_vf(j)%sf(i, q, l) + mom_3d(0, 0, k)%sf(i, q, l) - ! rhs_vf(j + 1)%sf(i, q, l) = & - ! rhs_vf(j + 1)%sf(i, q, l) + mom_3d(1, 0, k)%sf(i, q, l) - ! rhs_vf(j + 2)%sf(i, q, l) = & - ! rhs_vf(j + 2)%sf(i, q, l) + mom_3d(0, 1, k)%sf(i, q, l) - ! rhs_vf(j + 3)%sf(i, q, l) = & - ! rhs_vf(j + 3)%sf(i, q, l) + mom_3d(2, 0, k)%sf(i, q, l) - ! rhs_vf(j + 4)%sf(i, q, l) = & - ! rhs_vf(j + 4)%sf(i, q, l) + mom_3d(1, 1, k)%sf(i, q, l) - ! rhs_vf(j + 5)%sf(i, q, l) = & - ! rhs_vf(j + 5)%sf(i, q, l) + mom_3d(0, 2, k)%sf(i, q, l) - ! j = j + 6 - ! end do - - ! end do - ! end do - ! end do - ! else -!!$acc parallel loop collapse(3) gang vector default(present) - ! do l = 0, p - ! do k = 0, n - ! do j = 0, m - ! divu%sf(j, k, l) = 0d0 - ! divu%sf(j, k, l) = & - ! 5d-1/dx(j)*(q_prim_qp%vf(contxe + id)%sf(j + 1, k, l) - & - ! q_prim_qp%vf(contxe + id)%sf(j - 1, k, l)) - - ! end do - ! end do - ! end do - - ! ndirs = 1; if (n > 0) ndirs = 2; if (p > 0) ndirs = 3 - - ! if (id == ndirs) then - ! call s_compute_bubble_source(bub_adv_src, bub_r_src, bub_v_src, bub_p_src, bub_m_src, divu, nbub, & - ! q_cons_qp%vf(1:sys_size), q_prim_qp%vf(1:sys_size), t_step, id, rhs_vf) - ! end if - ! end if - !end if - - !if (monopole) then - ! ndirs = 1; if (n > 0) ndirs = 2; if (p > 0) ndirs = 3 - ! if (id == ndirs) then - ! call s_monopole_calculations(mono_mass_src, mono_mom_src, mono_e_src, & - ! q_cons_qp%vf(1:sys_size), q_prim_qp%vf(1:sys_size), t_step, id, & - ! rhs_vf) - ! end if - !end if - - !if (model_eqns == 3) then - ! !$acc parallel loop collapse(4) gang vector default(present) - ! do l = 0, p - ! do k = 0, n - ! do j = 0, m - ! do i = 1, num_fluids - ! rhs_vf(i + intxb - 1)%sf(j, k, l) = & - ! rhs_vf(i + intxb - 1)%sf(j, k, l) - 1d0/dx(j)* & - ! q_cons_qp%vf(i + advxb - 1)%sf(j, k, l)* & - ! q_prim_qp%vf(E_idx)%sf(j, k, l)* & - ! (flux_src_n(1)%vf(advxb)%sf(j, k, l) - & - ! flux_src_n(1)%vf(advxb)%sf(j - 1, k, l)) - ! end do - ! end do - ! end do - ! end do - !end if - - !if (any(Re_size > 0)) then - ! !$acc parallel loop collapse(3) gang vector default(present) - ! do l = 0, p - ! do k = 0, n - ! do j = 0, m - ! !$acc loop seq - ! do i = momxb, E_idx - ! rhs_vf(i)%sf(j, k, l) = & - ! rhs_vf(i)%sf(j, k, l) + 1d0/dx(j)* & - ! (flux_src_n(1)%vf(i)%sf(j - 1, k, l) & - ! - flux_src_n(1)%vf(i)%sf(j, k, l)) - ! end do - ! end do - ! end do - ! end do - !end if - elseif (id == 2) then ! RHS Contribution in y-direction =============================== ! Applying the Riemann fluxes @@ -1549,43 +1290,40 @@ contains if(barotropic) then !$acc parallel loop gang vector collapse(2) default(present) do k = -1, n+1 - do j = 0, m + do j = -1, m+1 flux_n(id)%vf(1)%sf(j,k,0) = q_prim_qp%vf(contxb)%sf(j,k,0) * q_prim_qp%vf(momxb+1)%sf(j,k,0) flux_n(id)%vf(2)%sf(j,k,0) = & q_prim_qp%vf(contxb)%sf(j,k,0) * & - (q_prim_qp%vf(momxb + 1)%sf(j,k,0)*q_prim_qp%vf(momxb)%sf(j,k,0)) + & - F_igr(2, j, k) + (q_prim_qp%vf(momxb + 1)%sf(j,k,0)*q_prim_qp%vf(momxb)%sf(j,k,0)) flux_n(id)%vf(3)%sf(j, k, 0) = q_prim_qp%vf(contxb)%sf(j,k,0) * & (q_prim_qp%vf(momxb + 1)%sf(j,k,0)**2d0) + & q_prim_qp%vf(e_idx)%sf(j,k,0) + & - F_igr(4, j, k) + F_igr(1, j, k) end do - end do - ! do i = 1, 3 - ! !$acc update host(flux_n(2)%vf(i)%sf) - ! end do - ! print *, "Flux rho Y", maxval(abs(flux_n(2)%vf(1)%sf( -1:m, -1:n+1, 0))) - ! print *, "Flux u Y", maxval(abs(flux_n(2)%vf(2)%sf( -1:m, -1:n+1, 0))) - ! print *, "Flux v Y", maxval(abs(flux_n(2)%vf(3)%sf( -1:m, -1:n+1, 0))) - - - ! !$acc update host(q_prim_qp%vf(momxb+1)%sf) - ! print *, "V MAX", maxval(abs(q_prim_qp%vf(momxb+1)%sf(0:m, -1:n+1, 0))) - !print *, "V LOC", maxloc(abs(q_prim_qp%vf(momxb+1)%sf(:, :, 0))) - - + end do !$acc parallel loop gang vector collapse(3) default(present) do i = 1, 3 do k = 0, n do j = 0, m - rhs_vf(i)%sf(j,k,0) = & - rhs_vf(i)%sf(j,k,0) + 1d0/(2d0*dy(k)) * & - ( flux_n(id)%vf(i)%sf(j,k-1,0) - & - flux_n(id)%vf(i)%sf(j,k+1,0) ) + if(lw == 1) then + rhs_vf(i)%sf(j,k,0) = rhs_vf(i)%sf(j,k,0) + & + 1d0/(2d0*dy(k)) * & + ( flux_n(id)%vf(i)%sf(j,k,0) - & + flux_n(id)%vf(i)%sf(j,k+1,0) + & + flux_n(id)%vf(i)%sf(j+1,k, 0) - & + flux_n(id)%vf(i)%sf(j+1,k+1, 0)) + else if(lw == 2) then + rhs_vf(i)%sf(j,k,0) = rhs_vf(i)%sf(j,k,0) + & + 1d0/(2d0*dy(k)) * & + ( flux_n(id)%vf(i)%sf(j-1,k-1,0) - & + flux_n(id)%vf(i)%sf(j-1,k,0) + & + flux_n(id)%vf(i)%sf(j,k-1, 0) - & + flux_n(id)%vf(i)%sf(j,k, 0)) + end if if(bcxb < -1) then @@ -1607,44 +1345,6 @@ contains end do end do end do - ! do i = 1, 3 - ! !$acc update host(rhs_vf(i)%sf) - ! end do - ! print *, "RHS RHO Y", maxval(abs(rhs_vf(1)%sf(0:m, 0:n, 0))) - ! print *, "RHS U Y", maxval(abs(rhs_vf(2)%sf(0:m, 0:n, 0))) - ! print *, "RHS V Y", maxval(abs(rhs_vf(3)%sf(0:m, 0:n, 0))) - - ! j = 149; k = 112 - ! print *, "VEL", q_prim_qp%vf(3)%sf(j, k, 0) - ! print *, "RHO", q_prim_qp%vf(1)%sf(j, k, 0) - ! print *, "PRES", q_prim_qp%vf(4)%sf(j, k, 0) - ! print *, "RHS IGR", rhs_igr(1:4, j, k) - ! print *, "JAC RHS SOL", jac_rhs_igr(2, j, k) , jac_igr(2, j, k) - ! print *, "F_IGR", F_igr(4, j, k) - - ! print *, "RHS ACTUAL", rhs_vf(3)%sf(j, k, 0) - ! print *, "FLUX ACTUAL", (1 / (2*dy(1)))* (flux_n(id)%vf(3)%sf(j, k-1, 0) - flux_n(id)%vf(3)%sf(j, k + 1, 0)) - ! print *, "FLUX_IGR_Y", (1 / (2*dy(1)))* (F_igr(4, j, k-1) - F_igr(4, j, k + 1)) - - ! k = 111 - ! print *, "PREV VEL", q_prim_qp%vf(3)%sf(j, k, 0) - ! print *, "PREV RHO", q_prim_qp%vf(1)%sf(j, k, 0) - ! print *, "PREV F", F_igr(4, j, k) - ! print *, "PREV PRES", q_prim_qp%vf(4)%sf(j, k, 0) - ! print *, "PREV JAC", jac_igr(2, j, k + 1), jac_igr(2, j, k) - ! print *, "PREV JAC RHS", jac_rhs_igr(2, j, k+1) , jac_rhs_igr(2, j, k) - ! print *, "PREV RHS", rhs_igr(4, j, k), rhs_igr(4, j, k-1), rhs_igr(4, j, k+1) - ! print *, "PREV DVDY", dvy_igr(j, k), dvy_igr(4, j, k-1) , dvy_igr(4, j, k+1) - - ! k = 113 - ! print *, "NEXT VEL", q_prim_qp%vf(3)%sf(j, k, 0) - ! print *, "NEXT RHO", q_prim_qp%vf(1)%sf(j, k, 0) - ! print *, "NEXT F", F_igr(4, j, k) - ! print *, "NEXT PRES", q_prim_qp%vf(4)%sf(j, k, 0) - ! print *, "NEXT JAC", jac_igr(2, j, k + 1), jac_igr(2, j, k) - ! print *, "NEXT JAC RHS", jac_rhs_igr(2, j, k+1) , jac_rhs_igr(2, j, k) - ! print *, "NEXT RHS", rhs_igr(4, j, k), rhs_igr(4, j, k-1), rhs_igr(4, j, k+1) - ! print *, "PREV DVDY", dvy_igr(j, k), dvy_igr(4, j, k-1) , dvy_igr(4, j, k+1) end if if(.not. barotropic) then !$acc parallel loop collapse(4) gang vector default(present) diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 43aa2dd9f..7987c81b4 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -346,9 +346,9 @@ contains do k = 0, n do j = 0, m q_cons_ts(2)%vf(i)%sf(j, k, l) = & - 0.25d0*(q_cons_ts(1)%vf(i)%sf(j - 1, k, l) + q_cons_ts(1)%vf(i)%sf(j + 1, k, l) + & - q_cons_ts(1)%vf(i)%sf(j , k + 1, l) + q_cons_ts(1)%vf(i)%sf(j ,k - 1, l)) & - + dt*rhs_vf(i)%sf(j, k, l) + 0.25d0*(q_cons_ts(1)%vf(i)%sf(j, k, l) + q_cons_ts(1)%vf(i)%sf(j + 1, k, l) + & + q_cons_ts(1)%vf(i)%sf(j , k + 1, l) + q_cons_ts(1)%vf(i)%sf(j + 1 ,k + 1, l)) & + + dt*0.5d0*rhs_vf(i)%sf(j, k, l) end do end do end do