From a6a16bb7d7711b3606f07539a4158483acc519c0 Mon Sep 17 00:00:00 2001 From: smaeyama Date: Wed, 3 Jul 2024 17:31:09 +0900 Subject: [PATCH 1/2] vp as indep. coordinate --- run/gkvp_namelist | 3 +- src/gkvp_advnc.f90 | 131 ++++++++++++++++++ src/gkvp_f0.56_bndry_tune_nec1.f90 | 209 +++++++++++++++++++++++++++++ src/gkvp_geom.f90 | 153 +++++++++++++++------ src/gkvp_header.f90 | 13 +- src/gkvp_out.f90 | 10 +- src/gkvp_set.f90 | 6 +- 7 files changed, 478 insertions(+), 47 deletions(-) diff --git a/run/gkvp_namelist b/run/gkvp_namelist index d75e94a..c314066 100644 --- a/run/gkvp_namelist +++ b/run/gkvp_namelist @@ -5,7 +5,8 @@ z_calc="cf4", art_diff=0.1d0, init_random=.false., - num_triad_diag=0, &end + num_triad_diag=0, + vp_coord=1, &end &triad mxt = 0, myt = 0/ &equib equib_type = "analytic", &end &run_n inum=%%%, diff --git a/src/gkvp_advnc.f90 b/src/gkvp_advnc.f90 index 3cd53a5..3640160 100644 --- a/src/gkvp_advnc.f90 +++ b/src/gkvp_advnc.f90 @@ -27,6 +27,10 @@ MODULE GKV_advnc use GKV_tips, only: tips_reality use GKV_geom, only: geom_increment_time +!> vp-mu + use GKV_bndry, only: bndry_shifts_m_e, bndry_shifts_m_f +!< vp-mu + implicit none private @@ -428,6 +432,12 @@ SUBROUTINE caldlt_linear( ff, psi, chi, dh ) deallocate( vb1o ) deallocate( vb2o ) +!> vp-mu + if ( vp_coord == 1 ) then + call literm_vp ( ff, psi, dh ) + end if +!< vp-mu + END SUBROUTINE caldlt_linear @@ -652,4 +662,125 @@ SUBROUTINE literm_zv ( ff, psi, im, lf ) END SUBROUTINE literm_zv +!> vp-mu +!-------------------------------------- + SUBROUTINE literm_vp ( ff, psi, lf ) +!-------------------------------------- +! z-derivative of 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), intent(in), & + dimension(-nx:nx,0:ny,-nz-nzb:nz-1+nzb,0:nm) :: psi + complex(kind=DP), intent(inout), & + dimension(-nx:nx,0:ny,-nz:nz-1,1:2*nv,0:nm) :: lf + + complex(kind=DP), dimension(:,:,:,:), allocatable :: psi_b + + real(kind=DP), dimension(-nz:nz-1) :: cefm + real(kind=DP) :: cs1, cs2 + integer :: mx, my, iz, iv, im + + + call clock_sta(1320) + ! call fapp_start("literm_perp",1320,1) + + allocate( psi_b(-nx:nx,0:ny,-nz-nzb:nz-1+nzb,0-nvb:nm+nvb) ) + + + cs1 = sgn(ranks) * Znum(ranks) / tau(ranks) + cs2 = sqrt( tau(ranks) / Anum(ranks) ) + + cefm(:) = 1._DP / ( 12._DP * dvp(:) ) * sqrt( tau(ranks) / Anum(ranks) ) + + + call bndry_shifts_m_e( psi, psi_b ) + call bndry_shifts_m_f( ff ) + + + if ( rankm /= 0 ) then + + + do im = 0, nm + do iv = 1, 2*nv + do iz = -nz, nz-1 + do my = ist_y, iend_y + do mx = -nx, nx + lf(mx,my,iz,iv,im) = lf(mx,my,iz,iv,im) & + - dpp(iz,iv,im) * cefm(iz) & + * ( - ff(mx,my,iz,iv,im+2) & + + 8._DP * ff(mx,my,iz,iv,im+1) & + - 8._DP * ff(mx,my,iz,iv,im-1) & + + ff(mx,my,iz,iv,im-2) & + + cs1 * fmx(iz,iv,im) * ( & + - psi_b(mx,my,iz,im+2) & + + 8._DP * psi_b(mx,my,iz,im+1) & + - 8._DP * psi_b(mx,my,iz,im-1) & + + psi_b(mx,my,iz,im-2) ) ) + end do + end do + end do + end do + end do + + else + + im = 1 + + do iv = 1, 2*nv + do iz = -nz, nz-1 + do my = ist_y, iend_y + do mx = -nx, nx + lf(mx,my,iz,iv,im) = lf(mx,my,iz,iv,im) & + - dpp(iz,iv,im) * cefm(iz) & + * ( - ff(mx,my,iz,iv,im+2) & + + 8._DP * ff(mx,my,iz,iv,im+1) & + - 8._DP * ff(mx,my,iz,iv,im-1) & + + ff(mx,my,iz,iv,im ) & + + cs1 * fmx(iz,iv,im) * ( & + - psi_b(mx,my,iz,im+2) & + + 8._DP * psi_b(mx,my,iz,im+1) & + - 8._DP * psi_b(mx,my,iz,im-1) & + + psi_b(mx,my,iz,im ) ) ) + end do + end do + end do + end do + + do im = 2, nm + do iv = 1, 2*nv + do iz = -nz, nz-1 + do my = ist_y, iend_y + do mx = -nx, nx + lf(mx,my,iz,iv,im) = lf(mx,my,iz,iv,im) & + - dpp(iz,iv,im) * cefm(iz) & + * ( - ff(mx,my,iz,iv,im+2) & + + 8._DP * ff(mx,my,iz,iv,im+1) & + - 8._DP * ff(mx,my,iz,iv,im-1) & + + ff(mx,my,iz,iv,im-2) & + + cs1 * fmx(iz,iv,im) * ( & + - psi_b(mx,my,iz,im+2) & + + 8._DP * psi_b(mx,my,iz,im+1) & + - 8._DP * psi_b(mx,my,iz,im-1) & + + psi_b(mx,my,iz,im-2) ) ) + end do + end do + end do + end do + end do + + + end if + + deallocate( psi_b ) + + ! call fapp_stop("literm_perp",1320,1) + call clock_end(1320) + + + END SUBROUTINE literm_vp +!< vp-mu + + + END MODULE GKV_advnc diff --git a/src/gkvp_f0.56_bndry_tune_nec1.f90 b/src/gkvp_f0.56_bndry_tune_nec1.f90 index ca0705f..b4fb825 100644 --- a/src/gkvp_f0.56_bndry_tune_nec1.f90 +++ b/src/gkvp_f0.56_bndry_tune_nec1.f90 @@ -23,6 +23,10 @@ MODULE GKV_bndry bndry_shifts_m_buffin, bndry_shifts_m_sendrecv, bndry_shifts_m_buffout, & bndry_vm_sendrecv_v2, bndry_zv_buffin_v2, bndry_zv_sendrecv_v2, bndry_zv_buffout_v2 +!> vp-mu + public bndry_shifts_m_e, bndry_shifts_m_f +!< vp-mu + CONTAINS @@ -1905,5 +1909,210 @@ SUBROUTINE bndry_zv_buffout_v2 ( zb2_bottom, zb2_top, vb2, ff ) END SUBROUTINE bndry_zv_buffout_v2 +!> 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_geom.f90 b/src/gkvp_geom.f90 index c246f10..f8ebd91 100644 --- a/src/gkvp_geom.f90 +++ b/src/gkvp_geom.f90 @@ -474,10 +474,27 @@ SUBROUTINE geom_init_kxkyzvm(lx, ly, eps_r_temp) dm = mmax / real( nprocm * ( nm+1 ) - 1, kind=DP ) ! --- equal spacing in vperp - do im = 0, nm - global_im = ( nm+1 ) * rankm + im - mu(im) = 0.5_DP * ( dm * real( global_im, kind=DP ) )**2 - end do +!> vp-mu + if( vp_coord == 1 ) then + + do im = 0, nm + global_im = ( nm+1 ) * rankm + im + do iz = -nz, nz-1 + vp(iz,im) = dm * real( global_im, kind=DP ) + end do + end do + + else + + do im = 0, nm + global_im = ( nm+1 ) * rankm + im + do iz = -nz, nz-1 + mu(iz,im) = 0.5_DP * ( dm * real( global_im, kind=DP ) )**2 + end do + end do + + end if +!< vp-mu do my = ist_y_g, iend_y_g @@ -1050,7 +1067,14 @@ SUBROUTINE geom_set_operators dpara(iz) = dz * q_0 * r_major do im = 0, nm - vp(iz,im) = sqrt( 2._DP * mu(im) )!* omg(iz) ) +!> vp-mu + if( vp_coord == 1 ) then + mu(iz,im) = 0.5_DP * vp(iz,im)**2 !/ omg(iz) + else + vp(iz,im) = sqrt( 2._DP * mu(iz,im) )!* omg(iz) ) + end if + dpp(iz,:,im) = 0._DP +!< vp-mu mir(iz,im) = 0._DP do iv = 1, 2*nv vdx(iz,iv,im) = 0._DP @@ -1058,7 +1082,7 @@ SUBROUTINE geom_set_operators vsy(iz,iv,im) = & - sgn(ranks) * tau(ranks) / Znum(ranks) & * ( R0_Ln(ranks) + R0_Lt(ranks) * ( 0.5_DP*vl(iv)**2 & - + omg(iz)*mu(im) - 1.5_DP ) ) + + omg(iz)*mu(iz,im) - 1.5_DP ) ) end do end do ! im loop ends @@ -1076,8 +1100,11 @@ SUBROUTINE geom_set_operators do im = 0, nm - vp(iz,im) = sqrt( 2._DP * mu(im) * omg(iz) ) - mir(iz,im) = mu(im) * eps_r / ( q_0 * r_major ) & +!> vp-mu + call geom_def_dpp( iz, im, domgdz ) +!< vp-mu + + mir(iz,im) = mu(iz,im) * eps_r / ( q_0 * r_major ) & * ( sin(zz_lab) & + eps_hor * lmmq * sin( lmmq * zz_lab - malpha ) & + eps_mor * lmmqm1 * sin( lmmqm1 * zz_lab - malpha ) & @@ -1085,7 +1112,7 @@ SUBROUTINE geom_set_operators do iv = 1, 2*nv vdx(iz,iv,im)= & - - ( vl(iv)**2 + omg(iz)*mu(im) ) * eps_rnew / r_major & + - ( vl(iv)**2 + omg(iz)*mu(iz,im) ) * eps_rnew / r_major & * ( 0._DP * ( rdeps00 + rdeps1_0 * cos( zz_lab ) & + rdeps2_10 * cos( lmmq * zz_lab - malpha ) & + rdeps1_10 * cos( lmmqm1 * zz_lab - malpha ) & @@ -1098,7 +1125,7 @@ SUBROUTINE geom_set_operators ) * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) vdy(iz,iv,im)= & - - ( vl(iv)**2 + omg(iz)*mu(im) ) * eps_rnew / r_major & + - ( vl(iv)**2 + omg(iz)*mu(iz,im) ) * eps_rnew / r_major & * ( 1._DP * ( rdeps00 + rdeps1_0 * cos( zz_lab ) & + rdeps2_10 * cos( lmmq * zz_lab - malpha ) & + rdeps1_10 * cos( lmmqm1 * zz_lab - malpha ) & @@ -1113,7 +1140,7 @@ SUBROUTINE geom_set_operators vsy(iz,iv,im) = & - sgn(ranks) * tau(ranks) / Znum(ranks) & * ( R0_Ln(ranks) + R0_Lt(ranks) * ( 0.5_DP*vl(iv)**2 & - + omg(iz)*mu(im) - 1.5_DP ) ) + + omg(iz)*mu(iz,im) - 1.5_DP ) ) end do @@ -1138,25 +1165,27 @@ SUBROUTINE geom_set_operators do im = 0, nm - vp(iz,im) = sqrt( 2._DP * mu(im) * omg(iz) ) +!> vp-mu + call geom_def_dpp( iz, im, domgdz ) +!< vp-mu - mir(iz,im) = mu(im) * (q_0/q_bar) * domgdz / ( omg(iz)*rootg(iz) ) + mir(iz,im) = mu(iz,im) * (q_0/q_bar) * domgdz / ( omg(iz)*rootg(iz) ) do iv = 1, 2*nv vdx(iz,iv,im) = & - ( vl(iv)**2 + omg(iz)*mu(im) ) / r_major & + ( vl(iv)**2 + omg(iz)*mu(iz,im) ) / r_major & * kkx & * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) vdy(iz,iv,im) = & - ( vl(iv)**2 + omg(iz)*mu(im) ) / r_major & + ( vl(iv)**2 + omg(iz)*mu(iz,im) ) / r_major & * kky & * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) vsy(iz,iv,im) = & - sgn(ranks) * tau(ranks) / Znum(ranks) & * ( R0_Ln(ranks) + R0_Lt(ranks) * ( 0.5_DP*vl(iv)**2 & - + omg(iz)*mu(im) - 1.5_DP ) ) & + + omg(iz)*mu(iz,im) - 1.5_DP ) ) & * (q_bar/q_0) end do @@ -1181,25 +1210,27 @@ SUBROUTINE geom_set_operators do im = 0, nm - vp(iz,im) = sqrt( 2._DP * mu(im) * omg(iz) ) +!> vp-mu + call geom_def_dpp( iz, im, domgdz ) +!< vp-mu - mir(iz,im) = mu(im) * domgdz / ( omg(iz)*rootg(iz) ) + mir(iz,im) = mu(iz,im) * domgdz / ( omg(iz)*rootg(iz) ) do iv = 1, 2*nv vdx(iz,iv,im)= & - ( vl(iv)**2 + omg(iz)*mu(im) ) / ( r_major*omg(iz) ) & + ( vl(iv)**2 + omg(iz)*mu(iz,im) ) / ( r_major*omg(iz) ) & * kkx & * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) vdy(iz,iv,im)= & - ( vl(iv)**2 + omg(iz)*mu(im) ) / ( r_major*omg(iz) ) & + ( vl(iv)**2 + omg(iz)*mu(iz,im) ) / ( r_major*omg(iz) ) & * kky & * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) vsy(iz,iv,im) = & - sgn(ranks) * tau(ranks) / Znum(ranks) & * ( R0_Ln(ranks) + R0_Lt(ranks) * ( 0.5_DP*vl(iv)**2 & - + omg(iz)*mu(im) - 1.5_DP ) ) + + omg(iz)*mu(iz,im) - 1.5_DP ) ) end do end do ! im loop ends @@ -1223,19 +1254,21 @@ SUBROUTINE geom_set_operators do im = 0, nm - vp(iz,im) = sqrt( 2._DP * mu(im) * omg(iz) ) +!> vp-mu + call geom_def_dpp( iz, im, domgdz ) +!< vp-mu - mir(iz,im) = mu(im) * domgdz / ( omg(iz)*rootg(iz) ) + mir(iz,im) = mu(iz,im) * domgdz / ( omg(iz)*rootg(iz) ) do iv = 1, 2*nv vdx(iz,iv,im) = & - ( vl(iv)**2 + omg(iz)*mu(im) ) / ( r_major*omg(iz) ) & + ( vl(iv)**2 + omg(iz)*mu(iz,im) ) / ( r_major*omg(iz) ) & * kkx & * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) vdy(iz,iv,im) = & - ( vl(iv)**2 + omg(iz)*mu(im) ) / ( r_major*omg(iz) ) & + ( vl(iv)**2 + omg(iz)*mu(iz,im) ) / ( r_major*omg(iz) ) & * kky & * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) & - real(ibprime,kind=DP) * vl(iv)**2 / r_major / omg(iz)**2 & ! grad-p (beta-prime) term @@ -1245,7 +1278,7 @@ SUBROUTINE geom_set_operators vsy(iz,iv,im) = & - sgn(ranks) * tau(ranks) / Znum(ranks) & * ( R0_Ln(ranks) + R0_Lt(ranks) * ( 0.5_DP*vl(iv)**2 & - + omg(iz)*mu(im) - 1.5_DP ) ) + + omg(iz)*mu(iz,im) - 1.5_DP ) ) end do end do ! im loop ends @@ -1268,18 +1301,20 @@ SUBROUTINE geom_set_operators do im = 0, nm - vp(iz,im) = sqrt( 2._DP * mu(im) * omg(iz) ) +!> vp-mu + call geom_def_dpp( iz, im, domgdz ) +!< vp-mu - mir(iz,im) = mu(im) * domgdz / ( omg(iz)*rootg(iz) ) + mir(iz,im) = mu(iz,im) * domgdz / ( omg(iz)*rootg(iz) ) do iv = 1, 2*nv vdx(iz,iv,im) = & - ( vl(iv)**2 + omg(iz)*mu(im) ) / ( r_major*omg(iz) ) & + ( vl(iv)**2 + omg(iz)*mu(iz,im) ) / ( r_major*omg(iz) ) & * kkx & * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) vdy(iz,iv,im) = & - ( vl(iv)**2 + omg(iz)*mu(im) ) / ( r_major*omg(iz) ) & + ( vl(iv)**2 + omg(iz)*mu(iz,im) ) / ( r_major*omg(iz) ) & * kky & * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) & - real(ibprime,kind=DP) * vl(iv)**2 / r_major / omg(iz)**2 & ! grad-p (beta-prime) term @@ -1289,7 +1324,7 @@ SUBROUTINE geom_set_operators vsy(iz,iv,im) = & - sgn(ranks) * tau(ranks) / Znum(ranks) & * ( R0_Ln(ranks) + R0_Lt(ranks) * ( 0.5_DP*vl(iv)**2 & - + omg(iz)*mu(im) - 1.5_DP ) ) + + omg(iz)*mu(iz,im) - 1.5_DP ) ) end do end do ! im loop ends @@ -1315,26 +1350,28 @@ SUBROUTINE geom_set_operators ! r_major = 1 is assumed as the equilibrium length unit ! B on the equatorial plane is also unity - vp(iz,im) = sqrt( 2._DP * mu(im) * omg(iz) ) +!> vp-mu + call geom_def_dpp( iz, im, domgdz ) +!< vp-mu - !mir(iz,im) = mu(im) * ub_dot_grdb - mir(iz,im) = mu(im) * domgdz / ( omg(iz)*rootg(iz) ) + !mir(iz,im) = mu(iz,im) * ub_dot_grdb + mir(iz,im) = mu(iz,im) * domgdz / ( omg(iz)*rootg(iz) ) do iv = 1, 2*nv vdx(iz,iv,im) = 0._DP !vdy(iz,iv,im) = & - ! ( vl(iv)**2 + omg(iz)*mu(im) ) & + ! ( vl(iv)**2 + omg(iz)*mu(iz,im) ) & ! * ( ub_crs_grdb / omg(iz)**2 ) * sqrt( gg(2,2) ) & ! * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) ! ion's vdy is negative y direction vdy(iz,iv,im) = & - ( vl(iv)**2 + omg(iz)*mu(im) ) / ( r_major*omg(iz) ) & + ( vl(iv)**2 + omg(iz)*mu(iz,im) ) / ( r_major*omg(iz) ) & * kky & * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) vsy(iz,iv,im)= & - sgn(ranks) * tau(ranks) / Znum(ranks) & * ( R0_Ln(ranks) + R0_Lt(ranks) * ( 0.5_DP*vl(iv)**2 & - + omg(iz)*mu(im) - 1.5_DP ) ) ! ion's vsy is negative y directuin + + omg(iz)*mu(iz,im) - 1.5_DP ) ) ! ion's vsy is negative y directuin end do end do ! im loop ends @@ -1360,7 +1397,7 @@ SUBROUTINE geom_set_operators do im = 0, nm do my = ist_y, iend_y do mx = -nx, nx - kmo = sqrt( 2._DP * ksq(mx,my,iz) * mu(im) / omg(iz) ) & + kmo = sqrt( 2._DP * ksq(mx,my,iz) * mu(iz,im) / omg(iz) ) & * dsqrt( tau(ranks)*Anum(ranks) ) / Znum(ranks) call math_j0( kmo, j0(mx,my,iz,im) ) call math_j1( kmo, j1(mx,my,iz,im) ) @@ -1392,7 +1429,14 @@ SUBROUTINE geom_set_operators if ( vel_rank == 0 ) then do iz = -nz, nz-1 - dvp(iz) = sqrt( 2._DP * (0.5_DP * dm**2) * omg(iz) ) +!> vp-mu +! dvp(iz) = sqrt( 2._DP * (0.5_DP * dm**2) * omg(iz) ) + if( vp_coord == 1 ) then + dvp(iz) = dm + else + dvp(iz) = sqrt( 2._DP * (0.5_DP * dm**2) * omg(iz) ) + end if +!< vp-mu end do end if call MPI_Bcast( dvp, 2*nz, MPI_DOUBLE_PRECISION, 0, & @@ -1401,7 +1445,7 @@ SUBROUTINE geom_set_operators do im = 0, nm do iv = 1, 2*nv do iz = -nz, nz-1 - fmx(iz,iv,im) = exp( - 0.5_DP * vl(iv)**2 - omg(iz) * mu(im) ) & + fmx(iz,iv,im) = exp( - 0.5_DP * vl(iv)**2 - omg(iz) * mu(iz,im) ) & / sqrt( twopi**3 ) end do end do @@ -1515,6 +1559,35 @@ SUBROUTINE geom_set_operators END SUBROUTINE geom_set_operators +!-------------------------------------- + SUBROUTINE geom_def_dpp( iz, im, domgdz ) +!-------------------------------------- + implicit none + + real(kind=DP), intent(in) :: domgdz + integer, intent(in) :: iz, im + + integer :: iv + + + if( vp_coord == 1 ) then + + mu(iz,im) = 0.5_DP * vp(iz,im)**2 / omg(iz) + do iv = 1, 2*nv + dpp(iz,iv,im) = 0.5_DP * vl(iv) * vp(iz,im) * domgdz & + / ( omg(iz) * omg(iz) * rootg(iz) ) + end do + + else + + vp(iz,im) = sqrt( 2._DP * mu(iz,im) * omg(iz) ) + dpp(iz,:,im) = 0._DP + + end if + + + END SUBROUTINE geom_def_dpp + !-------------------------------------- SUBROUTINE geom_reset_time(time_shearflow) !-------------------------------------- diff --git a/src/gkvp_header.f90 b/src/gkvp_header.f90 index f40f1db..6502bb0 100644 --- a/src/gkvp_header.f90 +++ b/src/gkvp_header.f90 @@ -45,7 +45,7 @@ MODULE GKV_header integer, parameter :: nxw = 20, nyw = 20 integer, parameter :: nx = 4, global_ny = 1 ! 2/3 de-aliasing rule - integer, parameter :: global_nz = 12, global_nv = 24, global_nm = 7 + integer, parameter :: global_nz = 12, global_nv = 24, global_nm = 31 integer, parameter :: nzb = 2, & ! the number of ghost grids in z nvb = 2 ! the number of ghost grids in v and m @@ -153,8 +153,12 @@ MODULE GKV_header real(kind=DP), dimension(0:ny) :: ky real(kind=DP), dimension(-nz:nz-1) :: zz, omg real(kind=DP), dimension(1:2*nv) :: vl - real(kind=DP), dimension(0:nm) :: mu - real(kind=DP), dimension(-nz:nz-1,0:nm) :: vp, mir +!> vp-mu +! real(kind=DP), dimension(0:nm) :: mu +! real(kind=DP), dimension(-nz:nz-1,0:nm) :: vp, mir + real(kind=DP), dimension(-nz:nz-1,0:nm) :: vp, mir, mu + real(kind=DP), dimension(-nz:nz-1,1:2*nv,0:nm) :: dpp ! the operator with 0.5*vl*vp*d/d(vp) +!< vp-mu real(kind=DP), dimension(-nz:nz-1) :: dvp real(kind=DP), dimension(-nz:nz-1) :: dpara, rootg @@ -195,6 +199,9 @@ MODULE GKV_header integer :: num_triad_diag +!> vp-mu + integer :: vp_coord +!< vp-mu !-------------------------------------- ! Parameters for numerical settings diff --git a/src/gkvp_out.f90 b/src/gkvp_out.f90 index 9da3197..812559e 100644 --- a/src/gkvp_out.f90 +++ b/src/gkvp_out.f90 @@ -1211,7 +1211,10 @@ SUBROUTINE write_moments ( ff, time ) do iz = -nz, nz-1 do my = ist_y, iend_y do mx = -nx, nx - wf(mx,my,iz,iv,im) = mu(im) * omg(iz) * ff(mx,my,iz,iv,im) * j0(mx,my,iz,im) +!> vp-mu +!!! wf(mx,my,iz,iv,im) = mu(im) * omg(iz) * ff(mx,my,iz,iv,im) * j0(mx,my,iz,im) + wf(mx,my,iz,iv,im) = mu(iz,im) * omg(iz) * ff(mx,my,iz,iv,im) * j0(mx,my,iz,im) +!< vp-mu end do end do end do @@ -1241,7 +1244,10 @@ SUBROUTINE write_moments ( ff, time ) do iz = -nz, nz-1 do my = ist_y, iend_y do mx = -nx, nx - wf(mx,my,iz,iv,im) = vl(iv) * mu(im) * omg(iz) * ff(mx,my,iz,iv,im) * j0(mx,my,iz,im) +!> vp-mu +!!! wf(mx,my,iz,iv,im) = vl(iv) * mu(im) * omg(iz) * ff(mx,my,iz,iv,im) * j0(mx,my,iz,im) + wf(mx,my,iz,iv,im) = vl(iv) * mu(iz,im) * omg(iz) * ff(mx,my,iz,iv,im) * j0(mx,my,iz,im) +!< vp-mu end do end do end do diff --git a/src/gkvp_set.f90 b/src/gkvp_set.f90 index c541732..9dc9b05 100644 --- a/src/gkvp_set.f90 +++ b/src/gkvp_set.f90 @@ -107,7 +107,10 @@ SUBROUTINE set_start namelist /cmemo/ memo namelist /calct/ calc_type, z_bound, z_filt, z_calc, art_diff, & - init_random, num_triad_diag +!> vp-mu +! init_random, num_triad_diag + init_random, num_triad_diag, vp_coord +!< vp-mu namelist /equib/ equib_type namelist /run_n/ inum, ch_res namelist /files/ f_log, f_hst, f_phi, f_fxv, f_cnt @@ -218,6 +221,7 @@ SUBROUTINE set_start write( olog, * ) "# Finite difference scheme in zz : ", trim(z_calc) write( olog, * ) "# Artificial diffusion in zz : ", art_diff write( olog, * ) "# Number of triad transfer diag. : ", num_triad_diag + write( olog, * ) "# Velocity space coordinate (=1 for vp, others for mu) : ", vp_coord write( olog, * ) "# Type of equib. : ", trim(equib_type) write( olog, * ) "" write( olog, * ) "# Run number = ", inum From 37b81b1b3a564c46d935cf48d2ccc631b51496bd Mon Sep 17 00:00:00 2001 From: smaeyama Date: Thu, 12 Dec 2024 16:46:09 +0900 Subject: [PATCH 2/2] 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)