From 37b81b1b3a564c46d935cf48d2ccc631b51496bd Mon Sep 17 00:00:00 2001 From: smaeyama Date: Thu, 12 Dec 2024 16:46:09 +0900 Subject: [PATCH] bndry, colliimp, dtc is updated for vp_coord. --- src/gkvp_bndry.f90 | 209 ++++++++++++++++++++++++++++++++++++++++++ src/gkvp_colliimp.f90 | 39 ++++++-- src/gkvp_dtc.f90 | 28 +++++- src/gkvp_geom.f90 | 2 + 4 files changed, 267 insertions(+), 11 deletions(-) diff --git a/src/gkvp_bndry.f90 b/src/gkvp_bndry.f90 index d3812b8..dee466c 100644 --- a/src/gkvp_bndry.f90 +++ b/src/gkvp_bndry.f90 @@ -25,6 +25,11 @@ MODULE GKV_bndry bndry_vm_buffin, bndry_vm_sendrecv, bndry_vm_buffout, & bndry_shifts_m_buffin, bndry_shifts_m_sendrecv, bndry_shifts_m_buffout +!> vp-mu + public bndry_shifts_m_e, bndry_shifts_m_f +!< vp-mu + + CONTAINS @@ -1490,5 +1495,209 @@ SUBROUTINE bndry_vm_buffout ( iz, vb2, mb2, ff ) END SUBROUTINE bndry_vm_buffout +!> vp-mu +!-------------------------------------- + SUBROUTINE bndry_shifts_m_e( ew, ew_b ) +!-------------------------------------- +! Shift communications in m direction + + complex(kind=DP), intent(in), & + dimension(-nx:nx,0:ny,-nz-nzb:nz-1+nzb,0:nm) :: ew + complex(kind=DP), intent(inout), & + dimension(-nx:nx,0:ny,-nz-nzb:nz-1+nzb,0-nvb:nm+nvb) :: ew_b + + complex(kind=DP), dimension(:,:,:,:), allocatable :: mb1, mb2 + + integer :: mx, my, iz, im + + allocate( mb1(-nx:nx,0:ny,-nz:nz-1,1:2*nvb) ) + allocate( mb2(-nx:nx,0:ny,-nz:nz-1,1:2*nvb) ) + + do im = 0, nm +!$OMP do schedule (dynamic) + do iz = -nz, nz-1 + do my = ist_y, iend_y + do mx = -nx, nx + ew_b(mx,my,iz,im) = ew(mx,my,iz,im) + end do + end do + end do +!$OMP end do nowait + end do + + call bndry_shifts_m_e_buffin( ew_b, mb1, mb2 ) + call bndry_shifts_m_e_sendrecv( mb1, mb2 ) + call bndry_shifts_m_e_buffout( mb2, ew_b ) + + deallocate( mb1 ) + deallocate( mb2 ) + + + END SUBROUTINE bndry_shifts_m_e + +!-------------------------------------- + SUBROUTINE bndry_shifts_m_f( ff ) +!-------------------------------------- + + complex(kind=DP), intent(inout), & + dimension(-nx:nx,0:ny,-nz-nzb:nz-1+nzb,1-nvb:2*nv+nvb,0-nvb:nm+nvb) :: ff + + complex(kind=DP), dimension(:,:,:,:,:), allocatable :: mb1, mb2 + + allocate( mb1(-nx:nx,0:ny,-nz:nz-1,1:2*nv,1:2*nvb) ) + allocate( mb2(-nx:nx,0:ny,-nz:nz-1,1:2*nv,1:2*nvb) ) + +!$OMP parallel default (none) & +!$OMP shared(ff,mb1,mb2) & + call bndry_shifts_m_buffin ( ff, mb1, mb2 ) +!$OMP barrier +!$OMP master + call bndry_shifts_m_sendrecv ( mb1, mb2 ) +!$OMP end master +!$OMP barrier + call bndry_shifts_m_buffout ( mb2, ff ) +!$OMP end parallel + + deallocate( mb1 ) + deallocate( mb2 ) + + END SUBROUTINE bndry_shifts_m_f + + +!-------------------------------------- + SUBROUTINE bndry_shifts_m_e_buffin( ff, mb1, mb2 ) +!-------------------------------------- +! Shift communications in m direction + + complex(kind=DP), intent(inout), & + dimension(-nx:nx,0:ny,-nz-nzb:nz-1+nzb,0-nvb:nm+nvb) :: ff + complex(kind=DP), intent(out), & + dimension(-nx:nx,0:ny,-nz:nz-1,1:2*nvb) :: mb1, mb2 + + integer :: mx, my, iz, im + + +!$OMP master + call clock_sta(1371) + ! call fapp_start("literm_shifts_bufferin",1371,1) +!$OMP end master + +! --- zero clear is required for rankv = 0, nprocv-1 and rankm = 0, nprocm-1 + do im = 1, 2*nvb +!$OMP do schedule (dynamic) + do iz = -nz, nz-1 + do my = ist_y, iend_y + do mx = -nx, nx + mb2(mx,my,iz,im) = ( 0._DP, 0._DP ) + end do + end do + end do +!$OMP end do nowait + end do + +!$OMP do schedule (dynamic) + do iz = -nz, nz-1 + do my = ist_y, iend_y + do mx = -nx, nx + do im = 1, nvb + mb1(mx,my,iz,im ) = ff(mx,my,iz, im-1) + mb1(mx,my,iz,im+nvb) = ff(mx,my,iz,nm-nvb+im) + end do + end do + end do + end do +!$OMP end do nowait + + +!$OMP master + ! call fapp_stop("literm_shifts_bufferin",1371,1) + call clock_end(1371) +!$OMP end master + + + END SUBROUTINE bndry_shifts_m_e_buffin + + +!-------------------------------------- + SUBROUTINE bndry_shifts_m_e_sendrecv( mb1, mb2 ) +!-------------------------------------- +! Shift communications in m direction + + complex(kind=DP), intent(in), & + dimension(-nx:nx,0:ny,-nz:nz-1,1:2*nvb) :: mb1 + complex(kind=DP), intent(inout), & + dimension(-nx:nx,0:ny,-nz:nz-1,1:2*nvb) :: mb2 + integer :: slngm + integer, dimension(4) :: ireq + integer, dimension(MPI_STATUS_SIZE,4) :: istatus + + + slngm = (2*nx+1)*(ny+1)*(2*nz)* nvb + + call clock_sta(1372) + ! call fapp_start("literm_shifts_sendrecv",1372,1) +! call MPI_sendrecv( mb1(-nx,0,-nz,1 ), slngm, MPI_DOUBLE_COMPLEX, imdn, 3, & +! mb2(-nx,0,-nz,nvb+1), slngm, MPI_DOUBLE_COMPLEX, imup, 3, & +! sub_comm_world, status, ierr_mpi ) +! +! call MPI_sendrecv( mb1(-nx,0,-nz,nvb+1), slngm, MPI_DOUBLE_COMPLEX, imup, 4, & +! mb2(-nx,0,-nz,1 ), slngm, MPI_DOUBLE_COMPLEX, imdn, 4, & +! sub_comm_world, status, ierr_mpi ) + + call MPI_irecv( mb2(-nx,0,-nz,nvb+1), slngm, MPI_DOUBLE_COMPLEX, imup, 3, & + sub_comm_world, ireq(1), ierr_mpi ) + call MPI_irecv( mb2(-nx,0,-nz, 1), slngm, MPI_DOUBLE_COMPLEX, imdn, 4, & + sub_comm_world, ireq(2), ierr_mpi ) + call MPI_isend( mb1(-nx,0,-nz, 1), slngm, MPI_DOUBLE_COMPLEX, imdn, 3, & + sub_comm_world, ireq(3), ierr_mpi ) + call MPI_isend( mb1(-nx,0,-nz,nvb+1), slngm, MPI_DOUBLE_COMPLEX, imup, 4, & + sub_comm_world, ireq(4), ierr_mpi ) + call MPI_waitall( 4, ireq, istatus, ierr_mpi ) + ! call fapp_stop("literm_shifts_sendrecv",1372,1) + call clock_end(1372) + + + END SUBROUTINE bndry_shifts_m_e_sendrecv + + +!-------------------------------------- + SUBROUTINE bndry_shifts_m_e_buffout( mb2, ff ) +!-------------------------------------- +! Shift communications in m direction + + complex(kind=DP), intent(in), & + dimension(-nx:nx,0:ny,-nz:nz-1,1:2*nvb) :: mb2 + complex(kind=DP), intent(inout), & + dimension(-nx:nx,0:ny,-nz-nzb:nz-1+nzb,0-nvb:nm+nvb) :: ff + + integer :: mx, my, iz, im + + +!$OMP master + call clock_sta(1373) + ! call fapp_start("literm_shifts_bufferout",1373,1) +!$OMP end master + +!$OMP do schedule (dynamic) + do iz = -nz, nz-1 + do my = ist_y, iend_y + do mx = -nx, nx + do im = 1, nvb + ff(mx,my,iz,-nvb-1+im) = mb2(mx,my,iz,im ) + ff(mx,my,iz,nm+im ) = mb2(mx,my,iz,im+nvb) + end do + end do + end do + end do +!$OMP end do nowait + +!$OMP master + ! call fapp_stop("literm_shifts_bufferout",1373,1) + call clock_end(1373) +!$OMP end master + + + END SUBROUTINE bndry_shifts_m_e_buffout +!< vp-mu END MODULE GKV_bndry diff --git a/src/gkvp_colliimp.f90 b/src/gkvp_colliimp.f90 index d2cbf5b..4a58ddf 100644 --- a/src/gkvp_colliimp.f90 +++ b/src/gkvp_colliimp.f90 @@ -35,7 +35,9 @@ MODULE GKV_colliimp ! end if !%%%%%%%%%%%% real(kind=DP), save :: gvl(1:2*global_nv) - real(kind=DP), save :: gmu(0:global_nm) +!> vp-mu + real(kind=DP), save :: gmu(0:global_nm,-nz:nz-1) +!< vp-mu real(kind=DP), save :: gvp(0:global_nm,-nz:nz-1) real(kind=DP), save :: gfmx(1:2*global_nv,0:global_nm,-nz:nz-1) real(kind=DP), save, & @@ -89,20 +91,35 @@ SUBROUTINE colliimp_set_param end do dm = vmax / real( nprocm * ( nm+1 ) - 1, kind=DP ) - do im = 0, global_nm - gmu(im) = 0.5_DP * ( dm * real( im, kind=DP ) )**2 - end do - do iz = -nz, nz-1 - do im = 0, global_nm - gvp(im,iz) = sqrt( 2._DP * gmu(im) * omg(iz) ) +!> vp-mu + if( vp_coord == 1 ) then + + do iz = -nz, nz-1 + do im = 0, global_nm + gvp(im,iz) = dm * real( im, kind=DP ) + gmu(im,iz) = 0.5_DP * gvp(im,iz)**2 / omg(iz) + end do end do - end do + + else + + do iz = -nz, nz-1 + do im = 0, global_nm + gmu(im,iz) = 0.5_DP * ( dm * real( im, kind=DP ) )**2 + gvp(im,iz) = sqrt( 2._DP * gmu(im,iz) * omg(iz) ) + end do + end do + + end if +!< vp-mu do iz = -nz, nz-1 do im = 0, global_nm do iv = 1, 2*global_nv - gfmx(iv,im,iz) = exp( - 0.5_DP * gvl(iv)**2 - omg(iz) * gmu(im) ) & +!> vp-mu + gfmx(iv,im,iz) = exp( - 0.5_DP * gvl(iv)**2 - omg(iz) * gmu(im,iz) ) & +!< vp-mu / sqrt( twopi**3 ) end do end do @@ -116,7 +133,9 @@ SUBROUTINE colliimp_set_param do iz = -nz, nz-1 do is = 0, ns-1 do im = 0, global_nm - kmo = sqrt( 2._DP * ksq(mx,my,iz) * gmu(im) / omg(iz) ) & +!> vp-mu + kmo = sqrt( 2._DP * ksq(mx,my,iz) * gmu(im,iz) / omg(iz) ) & +!< vp-mu * sqrt( tau(is)*Anum(is) ) / Znum(is) call math_j0( kmo, gj0(im,is,iz,ibuff) ) call math_j1( kmo, gj1(im,is,iz,ibuff) ) diff --git a/src/gkvp_dtc.f90 b/src/gkvp_dtc.f90 index 1aad21a..160bfef 100644 --- a/src/gkvp_dtc.f90 +++ b/src/gkvp_dtc.f90 @@ -48,6 +48,10 @@ SUBROUTINE dtc_init( lx, ly, vmax ) real(kind=DP) :: cs, dx, dy, kvd integer :: mx, my, iz, iv, im, is +!> vp-mu + real(kind=DP) :: dt_vp, dpp_max, dpp_max2 +!< vp-mu + ksq_max0 = 0._DP ksq_max = 0._DP @@ -104,6 +108,23 @@ SUBROUTINE dtc_init( lx, ly, vmax ) mir_max2 = max(mir_max2, 1.d-20) dt_vl = courant_num * dv / mir_max2 +!> vp-mu + dpp_max = 0._DP + do im = 0, nm + do iv = 1, 2*nv + do iz = -nz, nz-1 + if ( dpp_max < cs * abs(dpp(iz,iv,im)) / dvp(iz) ) then + dpp_max = cs * abs(dpp(iz,iv,im)) / dvp(iz) + end if + end do + end do + end do + call MPI_Allreduce( dpp_max, dpp_max2, 1, MPI_DOUBLE_PRECISION, & + MPI_MAX, MPI_COMM_WORLD, ierr_mpi ) + dpp_max2 = max(dpp_max2, 1.d-20) + dt_vp = courant_num / dpp_max2 * 0.25d0 ! a factor of 1/4 is added +!< vp-mu + if ( trim(col_type) == "LB" ) then nu_max = nu(ranks)*3._DP*dsqrt(pi)*ctauiv(ranks,ranks)/4._DP & @@ -180,7 +201,9 @@ SUBROUTINE dtc_init( lx, ly, vmax ) end if - dt_linear = min( dt_perp, dt_zz, dt_vl ) +!> vp-mu + dt_linear = min( dt_perp, dt_zz, dt_vl, dt_vp ) +!< vp-mu if (trim(time_advnc) == "rkg4") then @@ -233,6 +256,9 @@ SUBROUTINE dtc_init( lx, ly, vmax ) write( olog, * ) " # dt_perp = ", dt_perp write( olog, * ) " # dt_zz = ", dt_zz write( olog, * ) " # dt_vl = ", dt_vl +!> vp-mu + write( olog, * ) " # dt_vp = ", dt_vp +!< vp-mu write( olog, * ) " # dt_col = ", dt_col write( olog, * ) " # dt_linear = ", dt_linear write( olog, * ) " # dt_max = ", dt_max diff --git a/src/gkvp_geom.f90 b/src/gkvp_geom.f90 index f8ebd91..777d7d2 100644 --- a/src/gkvp_geom.f90 +++ b/src/gkvp_geom.f90 @@ -1559,6 +1559,7 @@ SUBROUTINE geom_set_operators END SUBROUTINE geom_set_operators +!> vp-mu !-------------------------------------- SUBROUTINE geom_def_dpp( iz, im, domgdz ) !-------------------------------------- @@ -1587,6 +1588,7 @@ SUBROUTINE geom_def_dpp( iz, im, domgdz ) END SUBROUTINE geom_def_dpp +!< vp-mu !-------------------------------------- SUBROUTINE geom_reset_time(time_shearflow)