diff --git a/README_for_namelist.txt b/README_for_namelist.txt index 470449b..9d93e2a 100644 --- a/README_for_namelist.txt +++ b/README_for_namelist.txt @@ -1,12 +1,13 @@ -NOTE for gkvp_f0.30 S. Maeyama March 2013 -Updated for gkvp_f0.40 M. Nakata June 2014 -Updated for gkvp_f0.45 M. Nakata July 2015 -Updated for gkvp_f0.46 S. Maeyama May 2016 -Updated for gkvp_f0.47 S. Maeyama Nov 2016 -Updated for gkvp_f0.48 S. Maeyama Dec 2016 -Updated for gkvp_f0.50 S. Maeyama Sep 2017 -Updated for gkvp_f0.55 M. Nakata Dec 2019 +Updated for gkvp_f0.62 S. Maeyama March 2023 Updated for gkvp_f0.61 S. Maeyama March 2021 +Updated for gkvp_f0.55 M. Nakata Dec 2019 +Updated for gkvp_f0.50 S. Maeyama Sep 2017 +Updated for gkvp_f0.48 S. Maeyama Dec 2016 +Updated for gkvp_f0.47 S. Maeyama Nov 2016 +Updated for gkvp_f0.46 S. Maeyama May 2016 +Updated for gkvp_f0.45 M. Nakata July 2015 +Updated for gkvp_f0.40 M. Nakata June 2014 +NOTE for gkvp_f0.30 S. Maeyama March 2013 %%% How to run the code %%% @@ -66,6 +67,7 @@ equib_type: "analytic" - Analytic helical field with the metrics in cylinder "vmec" - Tokamak/stellarator field from the VMEC code "eqdsk" - Tokamak field (MEUDAS/TOPICS or G-EQDSK) via IGS code "slab" - Shearless slab geometry + "ring" - Ring dipole geometry inum: current shot number @@ -156,6 +158,15 @@ del_c: mode connection phase in fluxtube model eps_r ~~ malpha : geometrical parameters such as safety factor, B-shear, etc. +&ring: parameters for ring dipole geometry + ! There is a ring current at R=a. The field line passing through (R,Z)=(R0,0) is picked up as a flux-tube domain. + ! The reference length is set to be R0 (not the ring current at R=a). + ! The reference magnetic field strength is B0 at (R,Z)=(R0,0). + +ring_a: = a / R0, which specify a flux tube of the ring dipole. + +kxmin: Minimum wavenumber in kx, valid only when equib_type == "ring" + &vmecp -- &bozxf : parameters for vmec equilibrium &igsp -- &igsf : parameters for tokamak (g-eqdsk) equilibrium diff --git a/Version_memo.txt b/Version_memo.txt index f6bde1b..6d11a21 100644 --- a/Version_memo.txt +++ b/Version_memo.txt @@ -1,3 +1,12 @@ +gkvp_f0.62 S. Maeyama Mar 2023 +1) equib_type = "ring" is added for ring dipole geometry. + +2) Rotating flux-tube model is implemented to treat equilibrium shearflows, + available for torus: equib_type = "s-alpha", "s-alpha-shift", "analytic", + "circMHD", "vmec", "igs". (But not available for "slab", "ring") + + + gkvp_f0.61 S. Maeyama Mar 2021 1) equib_type = "s-alpha-shift" is added. s-alpha model with Shafranov shift. @@ -7,7 +16,6 @@ gkvp_f0.61 S. Maeyama Mar 2021 - gkvp_f0.60 S. Maeyama Feb 2021 1) NetCDF4+parallel HDF5 is added for optional output of GKV. diff --git a/extra_tools/fig_stdout_20210316.tar.gz b/extra_tools/backup/fig_stdout_20210316.tar.gz similarity index 100% rename from extra_tools/fig_stdout_20210316.tar.gz rename to extra_tools/backup/fig_stdout_20210316.tar.gz diff --git a/extra_tools/fig_stdout_f0.62.tar.gz b/extra_tools/fig_stdout_f0.62.tar.gz new file mode 100644 index 0000000..48bf254 Binary files /dev/null and b/extra_tools/fig_stdout_f0.62.tar.gz differ diff --git a/lib/gkvp_math_portable.f90 b/lib/gkvp_math_portable.f90 index 5a8418e..a151358 100644 --- a/lib/gkvp_math_portable.f90 +++ b/lib/gkvp_math_portable.f90 @@ -5,6 +5,8 @@ MODULE GKV_math ! ! Update history of gkvp_set.f90 ! -------------- +! gkvp_f0.62 (S. Maeyama, Mar 2023) +! - Elliptic integrals math_eli1, math_eli2 are added. ! gkvp_f0.61 (S. Maeyama, Mar 2021) ! - random_seed is added for reproducibility. ! @@ -17,7 +19,9 @@ MODULE GKV_math private public math_j0, math_j0zero, math_j1, math_j2, & - math_i0, math_g0, math_random + math_i0, math_g0, math_random, & + math_eli1, math_eli2 + integer, parameter :: ifnc = 99 ! unit number preserved ! for this module @@ -158,6 +162,47 @@ SUBROUTINE math_g0( x, g0 ) END SUBROUTINE math_g0 + SUBROUTINE math_eli1( x, eli1 ) +! +! Complete elliptic integral of the first kind K +! + real(kind=DP), intent(in) :: x + real(kind=DP), intent(out) :: eli1 +! + real(kind=DP) :: sqrtx +! + sqrtx = sqrt(x) + if( 0._DP <= sqrtx .and. sqrtx < 1._DP ) then + eli1 = ellipk(sqrtx) + else + print *, "### math_eli1: x is out of range!" + end if + + return + + END SUBROUTINE math_eli1 + + SUBROUTINE math_eli2( x, eli2 ) +! +! 0th-order modified Bessel function +! + real(kind=DP), intent(in) :: x + real(kind=DP), intent(out) :: eli2 +! + real(kind=DP) :: sqrtx +! + sqrtx = sqrt(x) + if( 0._DP <= sqrtx .and. sqrtx < 1._DP ) then + eli2 = ellipe(sqrtx) + else + print *, "### math_eli2: x is out of range!" + end if + + return + + END SUBROUTINE math_eli2 + + SUBROUTINE math_random( rr ) ! @@ -825,6 +870,160 @@ function dbesi1(x) dbesi1 = y end function dbesi1 ! +FUNCTION ellipk(k) +! +! Complete elliptic integral of the 1st kind +! K(k) = \int_0^(pi/2) (1-k**2*sin(q)**2)**(-1/2) dq +! where 0<=k<1. +! +! Algorithm based on Jacobi's nome q expansion [Refs1-3] +! Implemented by S. Maeyama (July,2021) +! +! Ref1. T.Fukushima, "Fast computation of complete elliptic integrals and +! Jacobian elliptic functions", Celest. Mech. Dyn. Astr. 105, +! 305-328 (2009). DOI 10.1007/s10569-009-9228-z +! Ref2. J.Yamauchi, T.Uno, S.Ichimatsu, "Denshikeisanki notameno +! suuchikeisanhou III" (Baifukan, 1971), p.258. (Japanese) +! Ref3. http://www.totoha.net/fc2_mirror2/Complete_Elliptical_Integral.pdf +! + real(kind=8) :: ellipk + real(kind=8), intent(in) :: k + + real(kind=8), parameter :: halfpi=2.d0*atan(1.d0) + real(kind=8), parameter :: piinv=1.d0/(2.d0*halfpi) + real(kind=8) :: kp2,kp,sqrtkp,eps,eps4,q,q4,q9,q16,q25,t3 + + if (k<0.7d0) then + + kp2=1.d0-k*k + kp=sqrt(kp2) + sqrtkp=sqrt(kp) + eps=0.5d0*(1.d0-sqrtkp)/(1.d0+sqrtkp) + eps4=eps*eps*eps*eps + q=eps*(1.d0+eps4*(2.d0+eps4*(15.d0 & + +eps4*(150.d0+eps4*(1707.d0+eps4*(20910.d0+eps4*(268616.d0))))))) + q4=q*q*q*q + q9=q4*q4*q + q16=q4*q4*q4*q4 + q25=q16*q9 + t3=1.d0+2.d0*(q+q4+q9+q16+q25) + ellipk=halfpi*t3*t3 + + else + + kp2=k*k ! modify for large k + kp=sqrt(kp2) + sqrtkp=sqrt(kp) + eps=0.5d0*(1.d0-sqrtkp)/(1.d0+sqrtkp) + eps4=eps*eps*eps*eps + q=eps*(1.d0+eps4*(2.d0+eps4*(15.d0 & + +eps4*(150.d0+eps4*(1707.d0+eps4*(20910.d0+eps4*(268616.d0))))))) + q4=q*q*q*q + q9=q4*q4*q + q16=q4*q4*q4*q4 + q25=q16*q9 + t3=1.d0+2.d0*(q+q4+q9+q16+q25) + ellipk=halfpi*t3*t3 + ellipk=-piinv*ellipk*log(q) ! modify for large k + + end if + + return +END FUNCTION ellipk + + +FUNCTION ellipe(k) +! +! Complete elliptic integral of the 2nd kind +! E(k) = \int_0^(pi/2) (1-k**2*sin(q)**")**(1/2) dq +! where 0<=k<1. +! +! Algorithm based on Jacobi's nome q expansion [Refs1-3] +! Implemented by S. Maeyama (July,2021) +! +! Ref1. T.Fukushima, "Fast computation of complete elliptic integrals and +! Jacobian elliptic functions", Celest. Mech. Dyn. Astr. 105, +! 305-328 (2009). DOI 10.1007/s10569-009-9228-z +! Ref2. J.Yamauchi, T.Uno, S.Ichimatsu, "Denshikeisanki notameno +! suuchikeisanhou III" (Baifukan, 1971), p.258. (Japanese) +! Ref3. http://www.totoha.net/fc2_mirror2/Complete_Elliptical_Integral.pdf +! + real(kind=8) :: ellipe + real(kind=8), intent(in) :: k + + real(kind=8), parameter :: piover6=2.d0*atan(1.d0)/3.d0 + real(kind=8), parameter :: halfpi=2.d0*atan(1.d0) + real(kind=8), parameter :: piinv=1.d0/(2.d0*halfpi) + real(kind=8) :: kp2,kp,sqrtkp,eps,eps4,q,q2,q4,q6,q8,q9,& + q10,q12,q14,q16,q18,q20,q22,q25,t3,t3sq,fac,wq + real(kind=8) :: ellipkp, ellipk, ellipep + + if (k<0.7d0) then + + kp2=1.d0-k*k + kp=sqrt(kp2) + sqrtkp=sqrt(kp) + eps=0.5d0*(1.d0-sqrtkp)/(1.d0+sqrtkp) + eps4=eps*eps*eps*eps + q=eps*(1.d0+eps4*(2.d0+eps4*(15.d0 & + +eps4*(150.d0+eps4*(1707.d0+eps4*(20910.d0+eps4*(268616.d0))))))) + q2=q*q + q4=q2*q2 + q6=q4*q2 + q8=q4*q4 + q9=q8*q + q10=q4*q6 + q12=q6*q6 + q14=q6*q8 + q16=q8*q8 + q18=q9*q9 + q20=q10*q10 + q22=q20*q2 + q25=q16*q9 + t3=1.d0+2.d0*(q+q4+q9+q16+q25) !& + t3sq=t3*t3 + fac=piover6/t3sq + wq=q2+3.d0*q4+4.d0*q6+7.d0*q8+6.d0*q10+12.d0*q12+8.d0*q14 & + +15.d0*q16+13.d0*q18+18.d0*q20+12.d0*q22 + ellipe=fac*(1.d0+(1.d0+kp2)*t3sq*t3sq-24.d0*wq) + + else + + kp2=k*k ! modify for large k + kp=sqrt(kp2) + sqrtkp=sqrt(kp) + eps=0.5d0*(1.d0-sqrtkp)/(1.d0+sqrtkp) + eps4=eps*eps*eps*eps + q=eps*(1.d0+eps4*(2.d0+eps4*(15.d0 & + +eps4*(150.d0+eps4*(1707.d0+eps4*(20910.d0+eps4*(268616.d0))))))) + q2=q*q + q4=q2*q2 + q6=q4*q2 + q8=q4*q4 + q9=q8*q + q10=q4*q6 + q12=q6*q6 + q14=q6*q8 + q16=q8*q8 + q18=q9*q9 + q20=q10*q10 + q22=q20*q2 + q25=q16*q9 + t3=1.d0+2.d0*(q+q4+q9+q16+q25) + t3sq=t3*t3 + fac=piover6/t3sq + wq=q2+3.d0*q4+4.d0*q6+7.d0*q8+6.d0*q10+12.d0*q12+8.d0*q14 & + +15.d0*q16+13.d0*q18+18.d0*q20+12.d0*q22 + ellipep=fac*(1.d0+(1.d0+kp2)*t3sq*t3sq-24.d0*wq) ! modify for large k + ellipkp=halfpi*t3sq ! modify for large k + ellipk=-piinv*ellipkp*log(q) ! modify for large k + ellipe=ellipk+(halfpi-ellipep*ellipk)/ellipkp ! modify for large k + + end if + + return + +END FUNCTION ellipe END MODULE GKV_math diff --git a/run/Makefile b/run/Makefile index e28c901..71a0a2c 100644 --- a/run/Makefile +++ b/run/Makefile @@ -50,6 +50,7 @@ gkvp: $(SRC)gkvp_header.f90\ $(SRC)gkvp_tips.f90\ $(SRC)gkvp_vmecbzx.f90\ $(SRC)gkvp_igs.f90\ + $(SRC)gkvp_ring.f90\ $(SRC)gkvp_bndry.f90\ $(SRC)gkvp_colli.f90\ $(SRC)$(FFT).f90\ @@ -57,6 +58,7 @@ gkvp: $(SRC)gkvp_header.f90\ $(SRC)gkvp_colliimp.f90\ $(SRC)gkvp_freq.f90\ $(SRC)gkvp_zfilter.f90\ + $(SRC)gkvp_geom.f90\ $(SRC)gkvp_exb.f90\ $(SRC)gkvp_trans.f90\ $(SRC)gkvp_advnc.f90\ @@ -75,6 +77,7 @@ gkvp: $(SRC)gkvp_header.f90\ $(FC) $(FFLAGS) -c $(SRC)gkvp_tips.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_vmecbzx.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_igs.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_ring.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_bndry.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_colli.f90 $(FC) $(FFLAGS) -c $(SRC)$(FFT).f90 $(INC) @@ -82,6 +85,7 @@ gkvp: $(SRC)gkvp_header.f90\ $(FC) $(FFLAGS) -c $(SRC)gkvp_colliimp.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_freq.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_zfilter.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_geom.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_exb.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_trans.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_advnc.f90 @@ -101,6 +105,7 @@ gkvp: $(SRC)gkvp_header.f90\ gkvp_tips.o\ gkvp_vmecbzx.o\ gkvp_igs.o\ + gkvp_ring.o\ gkvp_bndry.o\ gkvp_colli.o\ $(FFT).o\ @@ -108,6 +113,7 @@ gkvp: $(SRC)gkvp_header.f90\ gkvp_colliimp.o\ gkvp_freq.o\ gkvp_zfilter.o\ + gkvp_geom.o\ gkvp_exb.o\ gkvp_trans.o\ gkvp_advnc.o\ diff --git a/run/backup/Makefile_flow b/run/backup/Makefile_flow index e28c901..71a0a2c 100644 --- a/run/backup/Makefile_flow +++ b/run/backup/Makefile_flow @@ -50,6 +50,7 @@ gkvp: $(SRC)gkvp_header.f90\ $(SRC)gkvp_tips.f90\ $(SRC)gkvp_vmecbzx.f90\ $(SRC)gkvp_igs.f90\ + $(SRC)gkvp_ring.f90\ $(SRC)gkvp_bndry.f90\ $(SRC)gkvp_colli.f90\ $(SRC)$(FFT).f90\ @@ -57,6 +58,7 @@ gkvp: $(SRC)gkvp_header.f90\ $(SRC)gkvp_colliimp.f90\ $(SRC)gkvp_freq.f90\ $(SRC)gkvp_zfilter.f90\ + $(SRC)gkvp_geom.f90\ $(SRC)gkvp_exb.f90\ $(SRC)gkvp_trans.f90\ $(SRC)gkvp_advnc.f90\ @@ -75,6 +77,7 @@ gkvp: $(SRC)gkvp_header.f90\ $(FC) $(FFLAGS) -c $(SRC)gkvp_tips.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_vmecbzx.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_igs.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_ring.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_bndry.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_colli.f90 $(FC) $(FFLAGS) -c $(SRC)$(FFT).f90 $(INC) @@ -82,6 +85,7 @@ gkvp: $(SRC)gkvp_header.f90\ $(FC) $(FFLAGS) -c $(SRC)gkvp_colliimp.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_freq.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_zfilter.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_geom.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_exb.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_trans.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_advnc.f90 @@ -101,6 +105,7 @@ gkvp: $(SRC)gkvp_header.f90\ gkvp_tips.o\ gkvp_vmecbzx.o\ gkvp_igs.o\ + gkvp_ring.o\ gkvp_bndry.o\ gkvp_colli.o\ $(FFT).o\ @@ -108,6 +113,7 @@ gkvp: $(SRC)gkvp_header.f90\ gkvp_colliimp.o\ gkvp_freq.o\ gkvp_zfilter.o\ + gkvp_geom.o\ gkvp_exb.o\ gkvp_trans.o\ gkvp_advnc.o\ diff --git a/run/backup/Makefile_fugaku b/run/backup/Makefile_fugaku index 8ae2a55..31a60ec 100644 --- a/run/backup/Makefile_fugaku +++ b/run/backup/Makefile_fugaku @@ -43,89 +43,99 @@ ifeq ($(FILEIO),gkvp_fileio_netcdf) LIB += -L$(NETCDF_DIR)/lib -lnetcdff -lnetcdf -lhdf5_hl -lhdf5 endif +OBJS = gkvp_header.o\ + gkvp_mpienv.o\ + $(MATH).o\ + gkvp_clock.o\ + $(FILEIO).o\ + gkvp_intgrl.o\ + gkvp_tips.o\ + gkvp_vmecbzx.o\ + gkvp_igs.o\ + gkvp_ring.o\ + gkvp_bndry.o\ + gkvp_colli.o\ + $(FFT).o\ + gkvp_fld.o\ + gkvp_colliimp.o\ + gkvp_freq.o\ + gkvp_zfilter.o\ + gkvp_geom.o\ + gkvp_exb.o\ + gkvp_trans.o\ + gkvp_advnc.o\ + gkvp_shearflow.o\ + gkvp_dtc.o\ + gkvp_out.o\ + gkvp_set.o\ + gkvp_main.o -gkvp: $(SRC)gkvp_header.f90\ - $(SRC)gkvp_mpienv.f90\ - $(MYL)$(MATH).f90\ - $(SRC)gkvp_clock.f90\ - $(SRC)$(FILEIO).f90\ - $(SRC)gkvp_intgrl.f90\ - $(SRC)gkvp_tips.f90\ - $(SRC)gkvp_vmecbzx.f90\ - $(SRC)gkvp_igs.f90\ - $(SRC)gkvp_bndry.f90\ - $(SRC)gkvp_colli.f90\ - $(SRC)$(FFT).f90\ - $(SRC)gkvp_fld.f90\ - $(SRC)gkvp_colliimp.f90\ - $(SRC)gkvp_freq.f90\ - $(SRC)gkvp_zfilter.f90\ - $(SRC)gkvp_exb.f90\ - $(SRC)gkvp_trans.f90\ - $(SRC)gkvp_advnc.f90\ - $(SRC)gkvp_shearflow.f90\ - $(SRC)gkvp_dtc.f90\ - $(SRC)gkvp_out.f90\ - $(SRC)gkvp_set.f90\ - $(SRC)gkvp_main.f90 +main: + (cp Makefile $(SRC); cd $(SRC); make gkvp) - $(FC) $(FFLAGS) -c $(SRC)gkvp_header.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_mpienv.f90 - $(FC) $(FFLAGS) -c $(MYL)$(MATH).f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_clock.f90 - $(FC) $(FFLAGS) -c $(SRC)$(FILEIO).f90 $(INC) - $(FC) $(FFLAGS) -c $(SRC)gkvp_intgrl.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_tips.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_vmecbzx.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_igs.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_bndry.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_colli.f90 - $(FC) $(FFLAGS) -c $(SRC)$(FFT).f90 $(INC) - $(FC) $(FFLAGS) -c $(SRC)gkvp_fld.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_colliimp.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_freq.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_zfilter.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_exb.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_trans.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_advnc.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_shearflow.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_dtc.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_out.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_set.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_main.f90 +gkvp: $(OBJS) + $(FC) $(FFLAGS) $(OBJS) -o $(PROG) $(LIB) + mv $(PROG) ../run/ - $(FC) $(FFLAGS) \ - gkvp_header.o\ - gkvp_mpienv.o\ - $(MATH).o\ - gkvp_clock.o\ - $(FILEIO).o\ - gkvp_intgrl.o\ - gkvp_tips.o\ - gkvp_vmecbzx.o\ - gkvp_igs.o\ - gkvp_bndry.o\ - gkvp_colli.o\ - $(FFT).o\ - gkvp_fld.o\ - gkvp_colliimp.o\ - gkvp_freq.o\ - gkvp_zfilter.o\ - gkvp_exb.o\ - gkvp_trans.o\ - gkvp_advnc.o\ - gkvp_shearflow.o\ - gkvp_dtc.o\ - gkvp_out.o\ - gkvp_set.o\ - gkvp_main.o\ - -o $(PROG) $(LIB) - - cp *.o *.mod *.$(OPTRPT) ../src/ - rm -f *.o *.mod *.$(OPTRPT) +#------------------------------> +gkvp_advnc.o : gkvp_advnc.f90 gkvp_geom.o gkvp_tips.o gkvp_zfilter.o gkvp_clock.o gkvp_bndry.o gkvp_colliimp.o gkvp_colli.o gkvp_exb.o gkvp_fld.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_bndry.o : gkvp_bndry.f90 gkvp_clock.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_clock.o : gkvp_clock.f90 gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_colli.o : gkvp_colli.f90 gkvp_bndry.o gkvp_clock.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_colliimp.o : gkvp_colliimp.f90 gkvp_fld.o $(MATH).o gkvp_clock.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_dtc.o : gkvp_dtc.f90 gkvp_colliimp.o gkvp_exb.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_exb.o : gkvp_exb.f90 gkvp_clock.o $(FFT).o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +$(FFT).o : $(FFT).f90 gkvp_clock.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< $(INC) +$(FILEIO).o : $(FILEIO).f90 gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< $(INC) +gkvp_fld.o : gkvp_fld.f90 gkvp_clock.o gkvp_intgrl.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_freq.o : gkvp_freq.f90 gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_geom.o : gkvp_geom.f90 gkvp_ring.o gkvp_igs.o gkvp_vmecbzx.o gkvp_intgrl.o $(MATH).o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_header.o : gkvp_header.f90 + $(FC) $(FFLAGS) -c $< +gkvp_igs.o : gkvp_igs.f90 gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_intgrl.o : gkvp_intgrl.f90 gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_main.o : gkvp_main.f90 gkvp_shearflow.o gkvp_tips.o gkvp_freq.o $(FFT).o gkvp_colliimp.o gkvp_advnc.o gkvp_fld.o gkvp_dtc.o gkvp_out.o gkvp_clock.o gkvp_set.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_mpienv.o : gkvp_mpienv.f90 gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_out.o : gkvp_out.f90 $(FILEIO).o gkvp_tips.o gkvp_dtc.o gkvp_colliimp.o gkvp_advnc.o gkvp_freq.o gkvp_trans.o gkvp_fld.o gkvp_intgrl.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_ring.o : gkvp_ring.f90 $(MATH).o + $(FC) $(FFLAGS) -c $< +gkvp_set.o : gkvp_set.f90 gkvp_geom.o $(FILEIO).o gkvp_tips.o gkvp_colliimp.o gkvp_colli.o gkvp_dtc.o gkvp_advnc.o gkvp_bndry.o gkvp_fld.o $(MATH).o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_shearflow.o : gkvp_shearflow.f90 gkvp_tips.o gkvp_fld.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_tips.o : gkvp_tips.f90 gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_trans.o : gkvp_trans.f90 $(FILEIO).o gkvp_exb.o gkvp_clock.o gkvp_intgrl.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_vmecbzx.o : gkvp_vmecbzx.f90 gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_vmecin.o : gkvp_vmecin.f90 gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_zfilter.o : gkvp_zfilter.f90 gkvp_clock.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +$(MATH).o : $(MYL)$(MATH).f90 $(MYL)Bessel0_Zeros.f90 gkvp_header.o + $(FC) $(FFLAGS) -c $< +#------------------------------< clean: - rm -f ../src/*.o ../src/*.mod ../src/*.$(OPTRPT) ./*.exe ./sub.q.*.o* \ + rm -f ../src/Makefile ../src/*.o ../src/*.mod ../src/*.$(OPTRPT) ./*.exe ./sub.q.*.o* \ ./*.o ./*.mod ./*.$(OPTRPT) ./*namelist.* ./sub.q.* clear: diff --git a/run/backup/Makefile_fugaku_old b/run/backup/Makefile_fugaku_old new file mode 100644 index 0000000..ec85633 --- /dev/null +++ b/run/backup/Makefile_fugaku_old @@ -0,0 +1,139 @@ +### Fujitsu Fortran Compiler ### +FC = mpifrtpx +FFLAGS = -Kfast,parallel # Optimization +FFLAGS += -X9 # Fortran95 +FFLAGS += -Koptmsg=2 -Nlst=t # Optimization report +FFLAGS += -fw # Suppress message +FFLAGS += -Kopenmp #-Nfjomplib # OpenMP +FFLAGS += -mcmodel=large # Static memory larger than 2GB +#FFLAGS += -Haefosux -NRtrap #-O0 # Debug +OPTRPT = 'lst' +#FFLAGS += -Nfjprof # Fujitsu profiler fapp +#FFLAGS += -Ksimd_nouse_multiple_structures # Specific option for compiler tcs1.2.26 to avoid slowing down GKV +#FFLAGS += -Knosch_pre_ra # Specific option for compiler tcs1.2.26 to avoid slowing down GKV + + +PROG = 'gkvp.exe' + +SRC = ../src/ +MYL = ../lib/ + +MATH = gkvp_math_portable + +FFT = gkvp_fft_fftw +### Usage of FFTW +ifeq ($(FFT),gkvp_fft_fftw) + ### FFTW-SVE + FFTW_DIR=/home/apps/r/OSS_CN/fftw-3.3.8/ + INC = -I$(FFTW_DIR)/include + LIB = -L$(FFTW_DIR)/lib64 -lfftw3 -lm -SSL2 + #### FFTW-SPACK (. /home/apps/oss/spack/share/spack/setup-env.sh; spack load fftw) ### + #FFTW_DIR=`spack location -i fftw` + #INC = -I$(FFTW_DIR)/include + #LIB = -L$(FFTW_DIR)/lib -lfftw3 -lm -SSL2 +endif + +FILEIO=gkvp_fileio_fortran +#FILEIO=gkvp_fileio_netcdf +### Usage of NetCDF (. /home/apps/oss/spack/share/spack/setup-env.sh; spack load netcdf-fortran%fj) ### +### Operation of NetCDF has not yet been checked on Fugaku, Jan 26 2021 +ifeq ($(FILEIO),gkvp_fileio_netcdf) + NETCDF_DIR=`spack location -i netcdf-fortran%fj` + INC += -I$(NETCDF_DIR)/include + LIB += -L$(NETCDF_DIR)/lib -lnetcdff -lnetcdf -lhdf5_hl -lhdf5 +endif + + +gkvp: $(SRC)gkvp_header.f90\ + $(SRC)gkvp_mpienv.f90\ + $(MYL)$(MATH).f90\ + $(SRC)gkvp_clock.f90\ + $(SRC)$(FILEIO).f90\ + $(SRC)gkvp_intgrl.f90\ + $(SRC)gkvp_tips.f90\ + $(SRC)gkvp_vmecbzx.f90\ + $(SRC)gkvp_igs.f90\ + $(SRC)gkvp_ring.f90\ + $(SRC)gkvp_bndry.f90\ + $(SRC)gkvp_colli.f90\ + $(SRC)$(FFT).f90\ + $(SRC)gkvp_fld.f90\ + $(SRC)gkvp_colliimp.f90\ + $(SRC)gkvp_freq.f90\ + $(SRC)gkvp_zfilter.f90\ + $(SRC)gkvp_geom.f90\ + $(SRC)gkvp_exb.f90\ + $(SRC)gkvp_trans.f90\ + $(SRC)gkvp_advnc.f90\ + $(SRC)gkvp_shearflow.f90\ + $(SRC)gkvp_dtc.f90\ + $(SRC)gkvp_out.f90\ + $(SRC)gkvp_set.f90\ + $(SRC)gkvp_main.f90 + + $(FC) $(FFLAGS) -c $(SRC)gkvp_header.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_mpienv.f90 + $(FC) $(FFLAGS) -c $(MYL)$(MATH).f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_clock.f90 + $(FC) $(FFLAGS) -c $(SRC)$(FILEIO).f90 $(INC) + $(FC) $(FFLAGS) -c $(SRC)gkvp_intgrl.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_tips.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_vmecbzx.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_igs.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_ring.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_bndry.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_colli.f90 + $(FC) $(FFLAGS) -c $(SRC)$(FFT).f90 $(INC) + $(FC) $(FFLAGS) -c $(SRC)gkvp_fld.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_colliimp.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_freq.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_zfilter.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_geom.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_exb.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_trans.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_advnc.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_shearflow.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_dtc.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_out.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_set.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_main.f90 + + $(FC) $(FFLAGS) \ + gkvp_header.o\ + gkvp_mpienv.o\ + $(MATH).o\ + gkvp_clock.o\ + $(FILEIO).o\ + gkvp_intgrl.o\ + gkvp_tips.o\ + gkvp_vmecbzx.o\ + gkvp_igs.o\ + gkvp_ring.o\ + gkvp_bndry.o\ + gkvp_colli.o\ + $(FFT).o\ + gkvp_fld.o\ + gkvp_colliimp.o\ + gkvp_freq.o\ + gkvp_zfilter.o\ + gkvp_geom.o\ + gkvp_exb.o\ + gkvp_trans.o\ + gkvp_advnc.o\ + gkvp_shearflow.o\ + gkvp_dtc.o\ + gkvp_out.o\ + gkvp_set.o\ + gkvp_main.o\ + -o $(PROG) $(LIB) + + cp *.o *.mod *.$(OPTRPT) ../src/ + rm -f *.o *.mod *.$(OPTRPT) + +clean: + rm -f ../src/*.o ../src/*.mod ../src/*.$(OPTRPT) ./*.exe ./sub.q.*.o* \ + ./*.o ./*.mod ./*.$(OPTRPT) ./*namelist.* ./sub.q.* + +clear: + rm -f ./*.o ./*.mod ./*.$(OPTRPT) ./*namelist.* ./sub.q.* + diff --git a/run/backup/Makefile_jfrs b/run/backup/Makefile_jfrs index 11ffec3..f771e28 100644 --- a/run/backup/Makefile_jfrs +++ b/run/backup/Makefile_jfrs @@ -61,6 +61,7 @@ gkvp: $(SRC)gkvp_header.f90\ $(SRC)gkvp_tips.f90\ $(SRC)gkvp_vmecbzx.f90\ $(SRC)gkvp_igs.f90\ + $(SRC)gkvp_ring.f90\ $(SRC)gkvp_bndry.f90\ $(SRC)gkvp_colli.f90\ $(SRC)$(FFT).f90\ @@ -68,6 +69,7 @@ gkvp: $(SRC)gkvp_header.f90\ $(SRC)gkvp_colliimp.f90\ $(SRC)gkvp_freq.f90\ $(SRC)gkvp_zfilter.f90\ + $(SRC)gkvp_geom.f90\ $(SRC)gkvp_exb.f90\ $(SRC)gkvp_trans.f90\ $(SRC)gkvp_advnc.f90\ @@ -87,6 +89,7 @@ gkvp: $(SRC)gkvp_header.f90\ $(FC) $(FFLAGS) -c $(SRC)gkvp_tips.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_vmecbzx.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_igs.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_ring.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_bndry.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_colli.f90 $(FC) $(FFLAGS) -c $(SRC)$(FFT).f90 $(INC) @@ -94,6 +97,7 @@ gkvp: $(SRC)gkvp_header.f90\ $(FC) $(FFLAGS) -c $(SRC)gkvp_colliimp.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_freq.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_zfilter.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_geom.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_exb.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_trans.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_advnc.f90 @@ -112,6 +116,7 @@ gkvp: $(SRC)gkvp_header.f90\ gkvp_tips.o\ gkvp_vmecbzx.o\ gkvp_igs.o\ + gkvp_ring.o\ gkvp_bndry.o\ gkvp_colli.o\ $(FFT).o\ @@ -119,6 +124,7 @@ gkvp: $(SRC)gkvp_header.f90\ gkvp_colliimp.o\ gkvp_freq.o\ gkvp_zfilter.o\ + gkvp_geom.o\ gkvp_exb.o\ gkvp_trans.o\ gkvp_advnc.o\ diff --git a/run/backup/Makefile_ps_sx b/run/backup/Makefile_ps_sx index 00307a5..d1a4262 100644 --- a/run/backup/Makefile_ps_sx +++ b/run/backup/Makefile_ps_sx @@ -47,6 +47,7 @@ gkvp: $(SRC)gkvp_header.f90\ $(SRC)gkvp_tips.f90\ $(SRC)gkvp_vmecbzx.f90\ $(SRC)gkvp_igs.f90\ + $(SRC)gkvp_ring.f90\ $(SRC)gkvp_f0.56_bndry_tune_nec1.f90\ $(SRC)gkvp_f0.56_colli_tune_nifs.f90\ $(SRC)$(FFT).f90\ @@ -56,8 +57,9 @@ gkvp: $(SRC)gkvp_header.f90\ $(SRC)gkvp_f0.56_zfilter_tune_nec1.f90\ $(SRC)gkvp_f0.56_exb_tune2r_0813.f90\ $(SRC)gkvp_trans.f90\ - $(SRC)gkvp_f0.56_advnc_tune_nec1.f90\ $(SRC)gkvp_shearflow.f90\ + $(SRC)gkvp_geom.f90\ + $(SRC)gkvp_advnc.f90\ $(SRC)gkvp_dtc.f90\ $(SRC)gkvp_out.f90\ $(SRC)gkvp_set.f90\ @@ -72,6 +74,7 @@ gkvp: $(SRC)gkvp_header.f90\ $(FC) $(FFLAGS) -c $(SRC)gkvp_tips.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_vmecbzx.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_igs.f90 + $(FC) $(FFRAGS) -c $(SRC)gkvp_ring.f90 $(FC) $(FFLAGS) $(FFLAGS_OMP3) -c $(SRC)gkvp_f0.56_bndry_tune_nec1.f90 $(FC) $(FFLAGS) $(FFLAGS_OMP2) -c $(SRC)gkvp_f0.56_colli_tune_nifs.f90 $(FC) $(FFLAGS) $(FFLAGS_OMP1) -c $(SRC)$(FFT).f90 $(INC) @@ -81,8 +84,9 @@ gkvp: $(SRC)gkvp_header.f90\ $(FC) $(FFLAGS) $(FFLAGS_OMP2) -c $(SRC)gkvp_f0.56_zfilter_tune_nec1.f90 $(FC) $(FFLAGS) $(FFLAGS_OMP1) -c $(SRC)gkvp_f0.56_exb_tune2r_0813.f90 $(FC) $(FFLAGS) $(FFLAGS_OMP1) -c $(SRC)gkvp_trans.f90 - $(FC) $(FFLAGS) $(FFLAGS_OMP3) -c $(SRC)gkvp_f0.56_advnc_tune_nec1.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_shearflow.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_geom.f90 + $(FC) $(FFLAGS) $(FFLAGS_OMP3) -c $(SRC)gkvp_advnc.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_dtc.f90 $(FC) $(FFLAGS) $(FFLAGS_OMP2) -c $(SRC)gkvp_out.f90 $(FC) $(FFLAGS) -c $(SRC)gkvp_set.f90 @@ -98,6 +102,7 @@ gkvp: $(SRC)gkvp_header.f90\ gkvp_tips.o\ gkvp_vmecbzx.o\ gkvp_igs.o\ + gkvp_ring.o\ gkvp_f0.56_bndry_tune_nec1.o\ gkvp_f0.56_colli_tune_nifs.o\ $(FFT).o\ @@ -107,8 +112,9 @@ gkvp: $(SRC)gkvp_header.f90\ gkvp_f0.56_zfilter_tune_nec1.o\ gkvp_f0.56_exb_tune2r_0813.o\ gkvp_trans.o\ - gkvp_f0.56_advnc_tune_nec1.o\ gkvp_shearflow.o\ + gkvp_geom.o\ + gkvp_advnc.o\ gkvp_dtc.o\ gkvp_out.o\ gkvp_set.o\ diff --git a/run/backup/Makefile_ubuntu b/run/backup/Makefile_ubuntu index 292007a..d39bc97 100644 --- a/run/backup/Makefile_ubuntu +++ b/run/backup/Makefile_ubuntu @@ -17,7 +17,8 @@ MATH = gkvp_math_portable FFT = gkvp_fft_fftw ifeq ($(FFT),gkvp_fft_fftw) - FFTW_DIR=`spack location -i fftw` + #FFTW_DIR=`spack location -i fftw` + FFTW_DIR=/usr INC += -I$(FFTW_DIR)/include LIB += -L$(FFTW_DIR)/lib -lfftw3 -lm endif @@ -30,90 +31,99 @@ ifeq ($(FILEIO),gkvp_fileio_netcdf) LIB += -L$(NETCDF_DIR)/lib -lnetcdff -lnetcdf -lhdf5_hl -lhdf5 endif +OBJS = gkvp_header.o\ + gkvp_mpienv.o\ + $(MATH).o\ + gkvp_clock.o\ + $(FILEIO).o\ + gkvp_intgrl.o\ + gkvp_tips.o\ + gkvp_vmecbzx.o\ + gkvp_igs.o\ + gkvp_ring.o\ + gkvp_bndry.o\ + gkvp_colli.o\ + $(FFT).o\ + gkvp_fld.o\ + gkvp_colliimp.o\ + gkvp_freq.o\ + gkvp_zfilter.o\ + gkvp_geom.o\ + gkvp_exb.o\ + gkvp_trans.o\ + gkvp_advnc.o\ + gkvp_shearflow.o\ + gkvp_dtc.o\ + gkvp_out.o\ + gkvp_set.o\ + gkvp_main.o -gkvp: $(SRC)gkvp_header.f90\ - $(SRC)gkvp_mpienv.f90\ - $(MYL)$(MATH).f90\ - $(SRC)gkvp_clock.f90\ - $(SRC)$(FILEIO).f90\ - $(SRC)gkvp_intgrl.f90\ - $(SRC)gkvp_tips.f90\ - $(SRC)gkvp_vmecbzx.f90\ - $(SRC)gkvp_igs.f90\ - $(SRC)gkvp_bndry.f90\ - $(SRC)gkvp_colli.f90\ - $(SRC)$(FFT).f90\ - $(SRC)gkvp_fld.f90\ - $(SRC)gkvp_colliimp.f90\ - $(SRC)gkvp_freq.f90\ - $(SRC)gkvp_zfilter.f90\ - $(SRC)gkvp_exb.f90\ - $(SRC)gkvp_trans.f90\ - $(SRC)gkvp_advnc.f90\ - $(SRC)gkvp_shearflow.f90\ - $(SRC)gkvp_dtc.f90\ - $(SRC)gkvp_out.f90\ - $(SRC)gkvp_set.f90\ - $(SRC)gkvp_main.f90 +main: + (cp Makefile $(SRC); cd $(SRC); make gkvp) - $(FC) $(FFLAGS) -c $(SRC)gkvp_header.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_mpienv.f90 - $(FC) $(FFLAGS) -c $(MYL)$(MATH).f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_clock.f90 - $(FC) $(FFLAGS) -c $(SRC)$(FILEIO).f90 $(INC) - $(FC) $(FFLAGS) -c $(SRC)gkvp_intgrl.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_tips.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_vmecbzx.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_igs.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_bndry.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_colli.f90 - $(FC) $(FFLAGS) -c $(SRC)$(FFT).f90 $(INC) - $(FC) $(FFLAGS) -c $(SRC)gkvp_fld.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_colliimp.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_freq.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_zfilter.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_exb.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_trans.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_advnc.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_shearflow.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_dtc.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_out.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_set.f90 - $(FC) $(FFLAGS) -c $(SRC)gkvp_main.f90 +gkvp: $(OBJS) + $(FC) $(FFLAGS) $(OBJS) -o $(PROG) $(LIB) + mv $(PROG) ../run/ - $(FC) $(FFLAGS) \ - gkvp_header.o\ - gkvp_mpienv.o\ - $(MATH).o\ - gkvp_clock.o\ - $(FILEIO).o\ - gkvp_intgrl.o\ - gkvp_tips.o\ - gkvp_vmecbzx.o\ - gkvp_igs.o\ - gkvp_bndry.o\ - gkvp_colli.o\ - $(FFT).o\ - gkvp_fld.o\ - gkvp_colliimp.o\ - gkvp_freq.o\ - gkvp_zfilter.o\ - gkvp_exb.o\ - gkvp_trans.o\ - gkvp_advnc.o\ - gkvp_shearflow.o\ - gkvp_dtc.o\ - gkvp_out.o\ - gkvp_set.o\ - gkvp_main.o\ - -o $(PROG) $(LIB) - -# cp *.o *.mod *.$(OPTRPT) ../src/ - cp *.o *.mod ../src/ - rm -f *.o *.mod *.$(OPTRPT) +#------------------------------> +gkvp_advnc.o : gkvp_advnc.f90 gkvp_geom.o gkvp_tips.o gkvp_zfilter.o gkvp_clock.o gkvp_bndry.o gkvp_colliimp.o gkvp_colli.o gkvp_exb.o gkvp_fld.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_bndry.o : gkvp_bndry.f90 gkvp_clock.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_clock.o : gkvp_clock.f90 gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_colli.o : gkvp_colli.f90 gkvp_bndry.o gkvp_clock.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_colliimp.o : gkvp_colliimp.f90 gkvp_fld.o $(MATH).o gkvp_clock.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_dtc.o : gkvp_dtc.f90 gkvp_colliimp.o gkvp_exb.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_exb.o : gkvp_exb.f90 gkvp_clock.o $(FFT).o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +$(FFT).o : $(FFT).f90 gkvp_clock.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< $(INC) +$(FILEIO).o : $(FILEIO).f90 gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< $(INC) +gkvp_fld.o : gkvp_fld.f90 gkvp_clock.o gkvp_intgrl.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_freq.o : gkvp_freq.f90 gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_geom.o : gkvp_geom.f90 gkvp_ring.o gkvp_igs.o gkvp_vmecbzx.o gkvp_intgrl.o $(MATH).o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_header.o : gkvp_header.f90 + $(FC) $(FFLAGS) -c $< +gkvp_igs.o : gkvp_igs.f90 gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_intgrl.o : gkvp_intgrl.f90 gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_main.o : gkvp_main.f90 gkvp_shearflow.o gkvp_tips.o gkvp_freq.o $(FFT).o gkvp_colliimp.o gkvp_advnc.o gkvp_fld.o gkvp_dtc.o gkvp_out.o gkvp_clock.o gkvp_set.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_mpienv.o : gkvp_mpienv.f90 gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_out.o : gkvp_out.f90 $(FILEIO).o gkvp_tips.o gkvp_dtc.o gkvp_colliimp.o gkvp_advnc.o gkvp_freq.o gkvp_trans.o gkvp_fld.o gkvp_intgrl.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_ring.o : gkvp_ring.f90 $(MATH).o + $(FC) $(FFLAGS) -c $< +gkvp_set.o : gkvp_set.f90 gkvp_geom.o $(FILEIO).o gkvp_tips.o gkvp_colliimp.o gkvp_colli.o gkvp_dtc.o gkvp_advnc.o gkvp_bndry.o gkvp_fld.o $(MATH).o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_shearflow.o : gkvp_shearflow.f90 gkvp_tips.o gkvp_fld.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_tips.o : gkvp_tips.f90 gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_trans.o : gkvp_trans.f90 $(FILEIO).o gkvp_exb.o gkvp_clock.o gkvp_intgrl.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_vmecbzx.o : gkvp_vmecbzx.f90 gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_vmecin.o : gkvp_vmecin.f90 gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +gkvp_zfilter.o : gkvp_zfilter.f90 gkvp_clock.o gkvp_mpienv.o gkvp_header.o + $(FC) $(FFLAGS) -c $< +$(MATH).o : $(MYL)$(MATH).f90 $(MYL)Bessel0_Zeros.f90 gkvp_header.o + $(FC) $(FFLAGS) -c $< +#------------------------------< clean: - rm -f ../src/*.o ../src/*.mod ../src/*.$(OPTRPT) ./*.exe ./sub.q.*.o* \ + rm -f ../src/Makefile ../src/*.o ../src/*.mod ../src/*.$(OPTRPT) ./*.exe ./sub.q.*.o* \ ./*.o ./*.mod ./*.$(OPTRPT) ./*namelist.* ./sub.q.* clear: diff --git a/run/backup/Makefile_ubuntu_old b/run/backup/Makefile_ubuntu_old new file mode 100644 index 0000000..4b81e37 --- /dev/null +++ b/run/backup/Makefile_ubuntu_old @@ -0,0 +1,128 @@ +#### gfortran + Open-MPI ### +FC = mpifort +FFLAGS = -mcmodel=medium -m64 -march=native -mtune=native -O3 -ffast-math # Optimization +#FFLAGS += -fopenmp # OpenMP +#FFLAGS += -Wall -Wextra -pedantic -fbacktrace \ +# -fbounds-check -Wuninitialized -ffpe-trap=invalid,zero,overflow # Debug +#INC = -I/usr/include +#LIB = -L/usr/lib + + +PROG = 'gkvp.exe' + +SRC = ../src/ +MYL = ../lib/ + +MATH = gkvp_math_portable + +FFT = gkvp_fft_fftw +ifeq ($(FFT),gkvp_fft_fftw) + #FFTW_DIR=`spack location -i fftw` + FFTW_DIR=/usr + INC += -I$(FFTW_DIR)/include + LIB += -L$(FFTW_DIR)/lib -lfftw3 -lm +endif + +FILEIO=gkvp_fileio_fortran +#FILEIO=gkvp_fileio_netcdf +ifeq ($(FILEIO),gkvp_fileio_netcdf) + NETCDF_DIR=`spack location -i netcdf-fortran` + INC += -I$(NETCDF_DIR)/include + LIB += -L$(NETCDF_DIR)/lib -lnetcdff -lnetcdf -lhdf5_hl -lhdf5 +endif + + +gkvp: $(SRC)gkvp_header.f90\ + $(SRC)gkvp_mpienv.f90\ + $(MYL)$(MATH).f90\ + $(SRC)gkvp_clock.f90\ + $(SRC)$(FILEIO).f90\ + $(SRC)gkvp_intgrl.f90\ + $(SRC)gkvp_tips.f90\ + $(SRC)gkvp_vmecbzx.f90\ + $(SRC)gkvp_igs.f90\ + $(SRC)gkvp_ring.f90\ + $(SRC)gkvp_bndry.f90\ + $(SRC)gkvp_colli.f90\ + $(SRC)$(FFT).f90\ + $(SRC)gkvp_fld.f90\ + $(SRC)gkvp_colliimp.f90\ + $(SRC)gkvp_freq.f90\ + $(SRC)gkvp_zfilter.f90\ + $(SRC)gkvp_geom.f90\ + $(SRC)gkvp_exb.f90\ + $(SRC)gkvp_trans.f90\ + $(SRC)gkvp_advnc.f90\ + $(SRC)gkvp_shearflow.f90\ + $(SRC)gkvp_dtc.f90\ + $(SRC)gkvp_out.f90\ + $(SRC)gkvp_set.f90\ + $(SRC)gkvp_main.f90 + + $(FC) $(FFLAGS) -c $(SRC)gkvp_header.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_mpienv.f90 + $(FC) $(FFLAGS) -c $(MYL)$(MATH).f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_clock.f90 + $(FC) $(FFLAGS) -c $(SRC)$(FILEIO).f90 $(INC) + $(FC) $(FFLAGS) -c $(SRC)gkvp_intgrl.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_tips.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_vmecbzx.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_igs.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_ring.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_bndry.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_colli.f90 + $(FC) $(FFLAGS) -c $(SRC)$(FFT).f90 $(INC) + $(FC) $(FFLAGS) -c $(SRC)gkvp_fld.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_colliimp.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_freq.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_zfilter.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_geom.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_exb.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_trans.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_advnc.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_shearflow.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_dtc.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_out.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_set.f90 + $(FC) $(FFLAGS) -c $(SRC)gkvp_main.f90 + + $(FC) $(FFLAGS) \ + gkvp_header.o\ + gkvp_mpienv.o\ + $(MATH).o\ + gkvp_clock.o\ + $(FILEIO).o\ + gkvp_intgrl.o\ + gkvp_tips.o\ + gkvp_vmecbzx.o\ + gkvp_igs.o\ + gkvp_ring.o\ + gkvp_bndry.o\ + gkvp_colli.o\ + $(FFT).o\ + gkvp_fld.o\ + gkvp_colliimp.o\ + gkvp_freq.o\ + gkvp_zfilter.o\ + gkvp_geom.o\ + gkvp_exb.o\ + gkvp_trans.o\ + gkvp_advnc.o\ + gkvp_shearflow.o\ + gkvp_dtc.o\ + gkvp_out.o\ + gkvp_set.o\ + gkvp_main.o\ + -o $(PROG) $(LIB) + +# cp *.o *.mod *.$(OPTRPT) ../src/ + cp *.o *.mod ../src/ + rm -f *.o *.mod *.$(OPTRPT) + +clean: + rm -f ../src/*.o ../src/*.mod ../src/*.$(OPTRPT) ./*.exe ./sub.q.*.o* \ + ./*.o ./*.mod ./*.$(OPTRPT) ./*namelist.* ./sub.q.* + +clear: + rm -f ./*.o ./*.mod ./*.$(OPTRPT) ./*namelist.* ./sub.q.* + diff --git a/run/backup/sub.q_ps_sx b/run/backup/sub.q_ps_sx index 46a59f8..ee004ca 100644 --- a/run/backup/sub.q_ps_sx +++ b/run/backup/sub.q_ps_sx @@ -53,10 +53,12 @@ MPI_procs=32 # number of MPI processes (= venode*8 for flat #PBS -v LANG=C +source /etc/profile.d/modules.sh + module load NECNLC-sx # module load NECNLC-mpi-sx ### For NetCDF -#module load netcdf-parallelIO-fortran-sx/4.5.2 +module load netcdf-parallelIO-fortran-sx ### Working directory diff --git a/run/gkvp_namelist b/run/gkvp_namelist index 76abf64..d75e94a 100644 --- a/run/gkvp_namelist +++ b/run/gkvp_namelist @@ -63,6 +63,9 @@ rdeps3_10= 0.d0, malpha = 0.d0, &end + &ring ring_a = 0.5d0, + kxmin = 0.05d0, &end + &vmecp s_input = 0.5d0, nss = 501, ntheta = 384, diff --git a/run/sub.q b/run/sub.q index ab115f4..6876642 100755 --- a/run/sub.q +++ b/run/sub.q @@ -68,7 +68,7 @@ DIR=%%DIR%% LDM=gkvp.exe NL=gkvp_namelist.%%% -#export XOS_MMM_L_PAGING_POLICY=demand:demand:demand # For Largepage +export XOS_MMM_L_PAGING_POLICY=demand:demand:demand # For Largepage export PLE_MPI_STD_EMPTYFILE="off" # Suppress stdout of filesize-0 diff --git a/src/gkvp_advnc.f90 b/src/gkvp_advnc.f90 index 4ccdb26..3cd53a5 100644 --- a/src/gkvp_advnc.f90 +++ b/src/gkvp_advnc.f90 @@ -5,6 +5,9 @@ MODULE GKV_advnc ! ! Update history of gkvp_advnc.f90 ! -------------- +! gkvp_f0.62 (S. Maeyama, Mar 2023) +! - Time-dependent metrics for rotating flux-tube model is implemented. +! See lines at "!%%% For shearflow rotating flux tube model %%%". ! gkvp_f0.57 (S. Maeyama, Oct 2020) ! - Version number f0.57 is removed from filename. ! - Unitialized access for padding iend_y metric_global_init + procedure :: xyz2rtq => metric_global_xyz2rtq + procedure :: rtq2xyz => metric_global_rtq2xyz + end type + + type metric_fourier + ! Metrics in flux coordinates at t=0, stored in Fourier coefficient + ! e.g., fourier_omg(kz) = \int omg(z) * exp(-i*kz*z) dz / \int dz + ! omg(z) = \sum_k fourier_omg(kz) * exp(i*kz*z) + ! Thus, omg(z) at arbitrary z is obtained by Fourier interpolation. + real(kind=DP), dimension(-global_nz:global_nz-1) :: kz + complex(kind=DP), dimension(-global_nz:global_nz-1) :: theta_tilde, omg + complex(kind=DP), dimension(-global_nz:global_nz-1) :: & + domgdr, domgdt, domgdq, grr, grt, grq, gtt, gtq, gqq, rootg_rtq + contains + procedure :: init => metric_fourier_init + procedure :: dft_rtq2coef => metric_fourier_dft_rtq2coef + end type + + type metric_local + ! Local metrics at any time t are stored. + ! They are updated with time integration, and used for solving + ! Gyrokinetic equation in the rotating flux tube model. + real(kind=DP), dimension(-nz:nz-1) :: zz ! The rotating flux tube coordinate (= z'') + real(kind=DP), dimension(-nz:nz-1) :: zz_labframe ! The flux-coordinate theta in the lab frame (= z''+t*gamma_e/s_hat) + real(kind=DP), dimension(-nz:nz-1) :: theta ! The geometrical poloidal angle theta_pol, not the flux-coordinate theta + real(kind=DP), dimension(-nz:nz-1) :: omg ! Magnetic field strength + real(kind=DP), dimension(-nz:nz-1) :: & + domgdx, domgdy, domgdz, gxx, gxy, gxz, gyy, gyz, gzz, rootg_xyz + real(kind=DP), dimension(-nz:nz-1) :: & + domgdr, domgdt, domgdq, grr, grt, grq, gtt, gtq, gqq, rootg_rtq + contains + procedure :: copy_global => metric_local_copy_global + procedure :: init => metric_local_init + procedure :: update => metric_local_update + procedure :: dft_coef2rtq => metric_local_dft_coef2rtq + procedure :: rtq2xyz => metric_local_rtq2xyz + end type + + type(metric_global), save :: mtr_global + type(metric_fourier), save :: mtr_fourier + type(metric_local), save :: mtr_local + + real(kind=DP), save :: cx, cy, cb + + +! for s-alpha model with Shafranov shift + real(kind=DP) :: p_total, dp_totaldx, beta_total, alpha_MHD + + real(kind=DP) :: r_major + + integer, parameter :: num_omtr = 13 +! real(kind=DP) :: metric_l(1:num_omtr,-nz:nz-1), metric_g(1:num_omtr,-global_nz:global_nz-1) + + real(kind=DP) :: s_hat + + real(kind=DP) :: eps_r + + real(kind=DP) :: lz, kxmin, kymin, dz, mmax, dm, del_c + real(kind=DP) :: z0, z0_l + integer :: n_tht, m_j + + real(kind=DP) :: rdeps00, eps_hor, lprd, mprd, lmmq, malpha + real(kind=DP) :: eps_mor, eps_por, lprdm1, lprdp1, lmmqm1, lmmqp1 + real(kind=DP) :: eps_rnew, rdeps1_0, rdeps1_10, rdeps2_10, rdeps3_10 + + real(kind=DP) :: s_input, s_0 ! radial label of fluxtube center + integer :: mc_type ! 0:Axisym., 1:Boozer, 2:Hamada + integer :: q_type ! 0:use q and s_hat value in confp, 1:calclated by IGS + integer :: isw, nss, ntheta, nzeta + real(kind=DP) :: phi_ax ! axisymmetric toroidal angle + +!sakano_ring-dipole st 202303 + real(kind=DP) :: ring_a +!sakano_ring-dipole end 202303 + + real(kind=DP) :: lz_l + + +CONTAINS + + +!-------------------------------------- + SUBROUTINE geom_read_nml +!-------------------------------------- + implicit none + + real(kind=DP) :: theta + real(kind=DP), dimension(0:ns-1) :: eta + real(kind=DP) :: domgdx, domgdy, domgdz + real(kind=DP), dimension(1:3,1:3) :: gg + integer :: iz, is, isw + + + namelist /physp/ R0_Ln, & ! R0/Lns + R0_Lt, & ! R0/Lts + nu, & ! factor for collision freq. in LB model + Anum, & ! mass number + Znum, & ! charge number + fcs, & ! charge-density fraction + sgn, & ! signs of charge + tau, & ! T-ratio Ts/T0, T0=reference ion temp. of ranks=1 + dns1, & ! initial perturbation amplitude + tau_ad, & ! Ti/Te for ITG-ae, Te/Ti for ETG-ai + lambda_i, & ! (Debye/rho_tp)^2 + beta, & ! mu0*ni*Ti/B^2 + ibprime, & ! flag for finite beta-prime effect on kvd + vmax, & ! maximum v_para in unit of v_ts + nx0 ! mode number for the initial perturbation + + namelist /rotat/ mach, uprime, gamma_e + + namelist /nperi/ n_tht, kymin, m_j, del_c + namelist /confp/ eps_r, eps_rnew, & + q_0, s_hat, & + lprd, mprd, eps_hor, eps_mor, eps_por, & + rdeps00, rdeps1_0, rdeps1_10, & + rdeps2_10, rdeps3_10, malpha +! namelist /vmecp/ q_0, rad_a, & +! R0_unit, r_edge, & +! b0b00, alpha_fix + namelist /vmecp/ s_input, nss, ntheta, nzeta + + namelist /igsp/ s_input, mc_type, q_type, nss, ntheta +!sakano_ring-dipole st 202303 + namelist /ring/ ring_a, kxmin +!sakano_ring-dipole end 202303 + + tau(:) = 1.0_DP + nu(:) = 0.002_DP + R0_Ln(:) = 2.5_DP + R0_Lt(:) = 7.5_DP + + + read(inml,nml=physp) + + + do is = 0, ns-1 + if( R0_Ln(is) /= 0._DP ) then + eta(is) = R0_Lt(is) / R0_Ln(is) + else + eta(is) = 1.d+20 + end if + end do + + + write( olog, * ) " # Physical parameters" + write( olog, * ) "" + write( olog, * ) " # r_major/L_ns = ", R0_Ln(:) + write( olog, * ) " # r_major/L_ts = ", R0_Lt(:) + write( olog, * ) " # eta = ", eta(:) + write( olog, * ) " # nu = ", nu(:) + write( olog, * ) " # A-number = ", Anum(:) + write( olog, * ) " # Z-number = ", Znum(:) + write( olog, * ) " # fcs = ", fcs(:) + write( olog, * ) " # sgn = ", sgn(:) + write( olog, * ) " # tau = ", tau(:) + write( olog, * ) " # dns1 = ", dns1(:) + write( olog, * ) " # tau_ad = ", tau_ad + write( olog, * ) " # lambda_i^2 = ", lambda_i + write( olog, * ) " # beta_i = ", beta + write( olog, * ) " # ibprime = ", ibprime + write( olog, * ) " # nx0 = ", nx0 + write( olog, * ) "" + + + mach = 0._DP + uprime = 0._DP + gamma_e = 0._DP + + read(inml,nml=rotat) + + write( olog, * ) " # Mean rotation parameters" + write( olog, * ) "" + write( olog, * ) " # Mach number = ", mach + write( olog, * ) " # uptime = ", uprime + write( olog, * ) " # gamma_ExB = ", gamma_e + write( olog, * ) "" + + + n_tht = 1 + + read(inml,nml=nperi) + + + if( trim(equib_type) == "slab") then + + read(inml,nml=confp) + + lprdm1 = 0._DP + lprdp1 = 0._DP + + lmmq = 0._DP + lmmqm1 = 0._DP + lmmqp1 = 0._DP + + q_0 = 1._DP ! For now, fixed q_0=1. Changing q_0 can extend parallel z-box size. + s_hat = 0._DP ! only shear less slab + eps_r = 1._DP + + eps_hor = 0._DP + lprd = 0._DP + mprd = 0._DP + malpha = 0._DP + + rdeps00 = 0._DP + eps_mor = 0._DP + eps_por = 0._DP + + write( olog, * ) " # Configuration parameters" + write( olog, * ) "" + write( olog, * ) " # q_0 = ", q_0 + write( olog, * ) " # s_hat = ", s_hat + write( olog, * ) " # eps_r = ", eps_r + write( olog, * ) "" + + write( olog, * ) " # eps_hor = ", eps_hor + write( olog, * ) " # lprd = ", lprd + write( olog, * ) " # mprd = ", mprd + write( olog, * ) " # malpha = ", malpha + write( olog, * ) " # rdeps00 = ", rdeps00 + + write( olog, * ) " # eps_mor = ", eps_mor + write( olog, * ) " # lprdm1 = ", lprdm1 + write( olog, * ) " # eps_por = ", eps_por + write( olog, * ) " # lprdp1 = ", lprdp1 + write( olog, * ) "" + + else if( trim(equib_type) == "analytic" .OR. & + trim(equib_type) == "s-alpha" .OR. & + trim(equib_type) == "s-alpha-shift" .OR. & + trim(equib_type) == "circ-MHD" ) then + + + read(inml,nml=confp) + + + lprdm1 = lprd - 1.0_DP + lprdp1 = lprd + 1.0_DP + + lmmq = lprd - mprd * q_0 + lmmqm1 = lprdm1 - mprd * q_0 + lmmqp1 = lprdp1 - mprd * q_0 + + + write( olog, * ) " # Configuration parameters" + write( olog, * ) "" + write( olog, * ) " # q_0 = ", q_0 + write( olog, * ) " # s_hat = ", s_hat + write( olog, * ) " # eps_r = ", eps_r + write( olog, * ) "" + + write( olog, * ) " # eps_hor = ", eps_hor + write( olog, * ) " # lprd = ", lprd + write( olog, * ) " # mprd = ", mprd + write( olog, * ) " # malpha = ", malpha + write( olog, * ) " # rdeps00 = ", rdeps00 + + write( olog, * ) " # eps_mor = ", eps_mor + write( olog, * ) " # lprdm1 = ", lprdm1 + write( olog, * ) " # eps_por = ", eps_por + write( olog, * ) " # lprdp1 = ", lprdp1 + write( olog, * ) "" + + + else if( trim(equib_type) == "vmec" ) then + + + read(inml,nml=confp) + + read(inml,nml=vmecp) + + call vmecbzx_boozx_read( nss, ntheta, nzeta ) + + isw = 0 + iz = 0 + call vmecbzx_boozx_coeff( isw, nss, ntheta, nzeta, s_input, iz, 0._DP, lz_l, & ! input + s_0, q_0, s_hat, eps_r, phi_ax, & ! output + omg(iz), rootg(iz), domgdx, domgdz, domgdy, & + gg(1,1), gg(1,2), gg(1,3), gg(2,2), & + gg(2,3), gg(3,3) ) + + write( olog, * ) " # Configuration parameters" + write( olog, * ) "" + write( olog, * ) " # r_major/L_ns = ", R0_Ln(:) + write( olog, * ) " # r_major/L_ts = ", R0_Lt(:) + write( olog, * ) " # eta = ", eta(:) + write( olog, * ) " # q_0 = ", q_0 + write( olog, * ) " # s_hat = ", s_hat + write( olog, * ) " # eps_r = ", eps_r + write( olog, * ) " # s_input, s_0 = ", s_input, s_0 + write( olog, * ) " # nss, ntheta, nzeta = ", nss, ntheta, nzeta + + + else if( trim(equib_type) == "eqdsk" ) then + + + read(inml,nml=confp) + + read(inml,nml=igsp) + + call igs_read( mc_type, nss, ntheta ) + + if ( q_type == 1 ) then + isw = 0 + iz = 0 + call igs_coeff( isw, mc_type, nss, ntheta, s_input, 0._DP, lz_l, & ! input + s_0, q_0, s_hat, eps_r, theta, & ! output + omg(iz), rootg(iz), domgdx, domgdz, domgdy, & + gg(1,1), gg(1,2), gg(1,3), gg(2,2), & + gg(2,3), gg(3,3) ) + end if + + write( olog, * ) " # Configuration parameters" + write( olog, * ) "" + write( olog, * ) " # r_major/L_ns = ", R0_Ln(:) + write( olog, * ) " # r_major/L_ts = ", R0_Lt(:) + write( olog, * ) " # eta = ", eta(:) + write( olog, * ) " # q_0 = ", q_0 + write( olog, * ) " # s_hat = ", s_hat + write( olog, * ) " # eps_r = ", eps_r + write( olog, * ) " # s_input, s_0 = ", s_input, s_0 + write( olog, * ) " # nss, ntheta = ", nss, ntheta + +!sakano_ring-dipole st 202303 + else if ( trim(equib_type) == "ring" ) then + + read(inml,nml=confp) + + read(inml,nml=ring) + + s_hat = 0._DP + + write( olog, * ) " # Configuration parameters for ring dipole configuration" + write( olog, * ) "" + write( olog, * ) " # s_hat = ", s_hat + write( olog, * ) " # kxmin = ", kxmin + write( olog, * ) " # ring_a = ", ring_a + write( olog, * ) " # eps_r = ", eps_r + write( olog, * ) " # q_0 = ", q_0 +!sakano_ring-dipole end 202303 + + else + + write( olog, * ) " # wrong choice of the equilibrium " + call flush(olog) + call MPI_Finalize(ierr_mpi) + stop + + end if + + END SUBROUTINE geom_read_nml + + +!-------------------------------------- + SUBROUTINE geom_init_kxkyzvm(lx, ly, eps_r_temp) +!-------------------------------------- + implicit none + real(kind=DP), intent(out) :: lx, ly, eps_r_temp + integer :: global_iv, global_im + integer :: mx, my, iz, iv, im + + eps_r_temp = eps_r + + if ( trim(equib_type) /= "ring" ) then + if (abs(s_hat) < 1.d-10) then ! When s_hat == ZERO + m_j = 0 + kxmin = kymin + else if (m_j == 0) then + kxmin = kymin + else + kxmin = abs(2._DP * pi * s_hat * kymin / real(m_j, kind=DP)) + end if + end if + lx = pi / kxmin + ly = pi / kymin + ! kymin=pi/ly=pi/[r_minor*pi/(q0*n_alp)]=q0*n_alp/r_minor + + lz = real( n_tht, kind=DP ) * pi ! total z-length + lz_l = lz / real( nprocz, kind=DP ) ! local z-length + + do mx = -nx, nx + kx(mx) = kxmin * real( mx, kind=DP ) + end do + + ky(:) = 0._DP + do my = ist_y_g, iend_y_g + ky(my-ist_y_g) = kymin * real( my, kind=DP ) + end do + + kxmin_g = kxmin + kymin_g = kymin + + z0 = - lz ! global lower boundary + z0_l = 2._DP * lz_l * real( rankz, kind=DP ) + z0 + ! local lower boundary + + dz = lz_l / real( nz, kind=DP ) + + do iz = -nz, nz-1 + zz(iz) = dz * real( iz + nz, kind=DP ) + z0_l + end do + + + dv = 2._DP * vmax / real( 2 * nv * nprocv -1, kind=DP ) + + do iv = 1, 2*nv + global_iv = 2 * nv * rankv + iv + vl(iv) = dv * ( real( global_iv - nv * nprocv - 1, kind=DP ) + 0.5_DP ) + end do + ! --- debug + ! write( olog, * ) " *** iv, vl " + ! do iv = 1, 2*nv + ! global_iv = 2 * nv * rankv + iv + ! write( olog, * ) iv, global_iv, vl(iv) + ! end do + ! write( olog, * ) "" + + mmax = vmax + 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 + + + do my = ist_y_g, iend_y_g + ck(my-ist_y_g) = exp( ui * 2._DP * pi * del_c & + * real( n_tht * my, kind=DP ) ) + dj(my-ist_y_g) = - m_j * n_tht * my + ! del_c = q_0*n_alp-int(q_0*n_alp) + ! m_j = 2*n_alp*q_d + end do + + + write( olog, * ) " # Numerical parameters" + write( olog, * ) "" + write( olog, * ) " # n_tht = ", n_tht + write( olog, * ) " # lx, ly, lz = ", lx, ly, lz + write( olog, * ) " # lz, z0 = ", lz, z0 + write( olog, * ) " # lz_l, z0_l = ", lz_l, z0_l + write( olog, * ) " # kxmin, kymin = ", kxmin, kymin + write( olog, * ) " # kxmax, kymax = ", kxmin*nx, kymin*global_ny + write( olog, * ) " # kperp_max = ", sqrt((kxmin*nx)**2+(kymin*global_ny)**2) + write( olog, * ) " # m_j, del_c = ", m_j, del_c + write( olog, * ) " # dz = ", dz + write( olog, * ) " # dv, vmax = ", dv, vmax + write( olog, * ) " # dm, mmax = ", dm, mmax + write( olog, * ) "" + + if (gamma_e == 0._DP) then + tlim_exb = 999999.d0 + else + tlim_exb = (kxmin*(nx-nx0))/(kymin*global_ny*abs(gamma_e)) + end if + write( olog, * ) " # ExB limit time tlim_exb = ", tlim_exb + write( olog, * ) " # for (mx=nx0,my=global_ny) initial perturbation: " + write( olog, * ) " # tlim_exb = kxmin*(nx-nx0)/(kymax*|gamma_e|)" + write( olog, * ) "" + + END SUBROUTINE geom_init_kxkyzvm + + +!-------------------------------------- + SUBROUTINE geom_init_metric +!-------------------------------------- + real(kind=DP) :: r_0 + real(kind=DP) :: wzz, theta, gomg + real(kind=DP) :: gdomgdx, gdomgdy, gdomgdz, & + ggxx, ggxy, ggxz, ggyy, ggyz, ggzz, grootg_xyz + real(kind=DP) :: gdomgdr, gdomgdt, gdomgdq, & + ggrr, ggrt, ggrq, ggtt, ggtq, ggqq, grootg_rtq + integer :: iz, is +!sakano_ring-dipole st 202303 + real(kind=DP) :: ub_dot_grdb, ub_crs_grdb +!sakano_ring-dipole end 202303 + + s_hat_g = s_hat + + !- zero clear - + gdomgdr = 0._DP + gdomgdt = 0._DP + gdomgdq = 0._DP + ggrr = 1._DP + ggrt = 0._DP + ggrq = 0._DP + ggtt = 1._DP + ggtq = 0._DP + ggqq = 1._DP + grootg_rtq = 1._DP + + do iz = -global_nz, global_nz-1 + + wzz = dz * iz + + if ( trim(equib_type) == "slab") then + !- Shearless slab geometry- + ! Consider translationally symmetric flux surface + ! (r,t,q)=(x_car,y_car,z_car). + ! GKV coordinates are + ! x = x_car, -lx<=x/gyroradius=B_ax, + ! where B_ax is the value at the magnetic axis. + ! cb = (psi_p(r))'/(cx*cy) = B_ax + ! Normalized = 1 and cb = 1 in the B_ax unit. + !- + r_major = 1._DP ! Major radius of magnetic axis in the R0 unit + r_0 = r_major * eps_r ! Minor radius of flux-tube center + cx = 1._DP + cy = r_0/q_0 + cb = 1._DP + + q_bar = q_0 + theta = wzz + gomg = 1._DP & + - eps_r * ( cos( wzz ) & + + eps_hor * cos( lmmq * wzz - malpha ) & + + eps_mor * cos( lmmqm1 * wzz - malpha ) & + + eps_por * cos( lmmqp1 * wzz - malpha ) ) + !- Metrics in GKV coordinates (x,y,z) + gdomgdz = eps_r * ( sin(wzz) & + + eps_hor * lmmq * sin( lmmq * wzz - malpha ) & + + eps_mor * lmmqm1 * sin( lmmqm1 * wzz - malpha ) & + + eps_por * lmmqp1 * sin( lmmqp1 * wzz - malpha ) ) + gdomgdy = - eps_rnew / r_major * ( & + - ( sin( wzz ) & + + eps_hor * lprd * sin( lmmq * wzz - malpha ) & + + eps_mor * lprdm1 * sin( lmmqm1 * wzz - malpha ) & + + eps_por * lprdp1 * sin( lmmqp1 * wzz - malpha ) & + ) - (-1._DP/eps_r) * gdomgdz ) + gdomgdx = eps_rnew / r_major * ( & + - ( & + rdeps00 & + + rdeps1_0 * cos( wzz ) & + + rdeps2_10 * cos( lmmq * wzz - malpha ) & + + rdeps1_10 * cos( lmmqm1 * wzz - malpha ) & + + rdeps3_10 * cos( lmmqp1 * wzz - malpha ) & + + s_hat * wzz * ( sin( wzz ) & + + eps_hor * lprd * sin( lmmq * wzz - malpha ) & + + eps_mor * lprdm1 * sin( lmmqm1 * wzz - malpha ) & + + eps_por * lprdp1 * sin( lmmqp1 * wzz - malpha ) ) & + ) - (-s_hat*wzz/eps_r) * gdomgdz ) + ggxx = 1._DP + ggxy = s_hat*wzz + ggxz = 0._DP + ggyy = 1._DP + (s_hat*wzz)**2 + ggyz = 1._DP/r_0 + ggzz = 1._DP/r_0**2 + grootg_xyz = q_0*r_major/gomg + !- Metrics in flux coordinates (r,theta,zeta) + !ggrr = 1._DP + !ggrt = 0._DP + !ggrq = 0._DP + !ggtt = 1._DP/r_0**2 + !ggtq = 0._DP + !ggqq = 0._DP ! /=1._DP/r_major**2, because of large-aspect ratio approximation + !grootg_rtq = r_0*r_major/gomg + + + else if( trim(equib_type) == "s-alpha" .or. trim(equib_type) == "s-alpha-shift" ) then + !- Analytic model of large-aspect-ratio tokamak system - + ! Consider concentric circular flux surface (r,theta,zeta). + ! GKV coordinates are + ! x = cx*(r-r0) + ! y = cy*(q(r)*theta-zeta) + ! z = theta + ! with cx=1, cy=cx*r0/q0. + ! In the large-aspect ratio limit, the geometrical length + ! in the field-aligned direction is dpara=q*r_major*dz. + ! r_major = 1 in the R0 unit. + ! Finite aspect ratio eps_r = r0/R0 is retained only in + ! magnetic field omg, domgdx, domgdy, domgdz, but not for metrics. + ! Flux-surface averaged magnetic field is =B_ax, + ! where B_ax is the value at the magnetic axis. + ! cb = (psi_p(r))'/(cx*cy) = B_ax + ! Normalized = 1 and cb = 1 in the B_ax unit. + !- + r_major = 1._DP ! Major radius of magnetic axis in the R0 unit + r_0 = r_major * eps_r ! Minor radius of flux-tube center + cx = 1._DP + cy = r_0/q_0 + cb = 1._DP + + if (trim(equib_type) == "s-alpha") then + !--- s-alpha model without Shafranov shift - + alpha_MHD = 0._DP + else if (trim(equib_type) == "s-alpha-shift") then + !--- s-alpha model with Shafranov shift ---- + p_total = 0._DP + dp_totaldx = 0._DP + beta_total = 0._DP + do is = 0, ns-1 + p_total = p_total + fcs(is) * tau(is) / Znum(is) + dp_totaldx = dp_totaldx - fcs(is) * tau(is) / Znum(is) * (R0_Ln(is) + R0_Lt(is)) + beta_total = beta_total + 2._DP * beta * fcs(is) * tau(is) / Znum(is) + end do + alpha_MHD = - q_0**2 * r_major * beta_total * dp_totaldx / p_total + end if + q_bar = q_0 + theta = wzz + gomg = 1._DP - eps_r * cos( theta ) ! s-alpha with eps-expansion + !!!!gomg = 1._DP / (1._DP + eps_r * cos( theta )) ! for benchmark + !- Metrics in GKV coordinates (x,y,z) + gdomgdz = eps_r * sin( theta ) + !!!!gdomgdz = eps_r * sin( theta ) * omg(iz)**2 ! for benchmark + gdomgdx = - cos( theta ) / r_major + gdomgdy = 0._DP + ggxx = 1._DP + ggxy = s_hat*wzz - alpha_MHD*sin(wzz) ! with Shafranov shift + ggxz = 0._DP + ggyy = 1._DP + (s_hat*wzz - alpha_MHD*sin(wzz))**2 ! with Shafranov shift + ggyz = 1._DP/r_0 + ggzz = 1._DP/r_0**2 + grootg_xyz = q_0*r_major/gomg + !- Metrics in flux coordinates (r,theta,zeta) + !gdomgdr = - cos( theta ) / r_major + !gdomgdt = eps_r * sin( theta ) + !gdomgdq = 0._DP + !ggrr = 1._DP + !ggrt = 0._DP + !ggrq = 0._DP + sin(wzz)*alpha_MHD*q_0/r_0 ! with Shafranov shift + !ggtt = 1._DP/r_0**2 + !ggtq = 0._DP + !ggqq = 0._DP & ! /=1._DP/r_major**2, because of large-aspect ratio approximation + ! + sin(wzz)*(alpha_MHD*q_0/r_0)**2 ! with Shafranov shift + !grootg_rtq = r_0*r_major/gomg + + + else if( trim(equib_type) == "circ-MHD" ) then + !- Circular MHD equilibrium - + ! [Ref.] X. Lapillonne, et al., Phys. Plasmas 16, 032308 (2009). + ! + ! Consider concentric circular flux surface (r,theta,zeta). + ! GKV coordinates are + ! x = cx*(r-r0) + ! y = cy*(q(r)*theta-zeta) + ! z = theta + ! with cx=1, cy=cx*r0/q0. + ! In the large-aspect ratio limit, the geometrical length + ! in the field-aligned direction is dpara=q*r_major*dz. + ! r_major = 1 in the R0 unit. + ! In contrast to the s-alpha model, finite aspect ratio eps_r = r0/R0 + ! is retained in both of magnetic field and metrics. + ! Difference between the flux-surface averaged magnetic field + ! and the value at the magnetic axis B_ax also appears. + ! cb = (psi_p(r))'/(cx*cy) = B_ax + ! Normalized omg = B(z)/B_ax and cb = 1 in the B_ax unit. + !- + r_major = 1._DP ! in the R0 unit + r_0 = r_major * eps_r ! Minor radius of flux-tube center + cx = 1._DP + cy = r_0/q_0 + cb = 1._DP + + theta = 2._DP*atan( sqrt( (1._DP+eps_r)/(1._DP-eps_r) ) & + * tan(wzz/2._DP) ) + q_bar = dsqrt( 1._DP - eps_r**2 )*q_0 + gomg = sqrt( q_bar**2 + eps_r**2 ) & + / ( 1._DP + eps_r*cos( theta ) ) / q_bar + !- Metrics in GKV coordinates (x,y,z) + gdomgdz = eps_r * sin(theta) * sqrt( q_bar**2 + eps_r**2 ) & + / ( 1._DP + eps_r * cos( theta ) )**2 & + / ( 1._DP - eps_r * cos( wzz) ) / q_0 + gdomgdx = -( cos(theta) & + - eps_r*(1._DP-s_hat+eps_r**2*q_0**2/q_bar**2) & + *(1._DP+eps_r*cos(theta))/(q_bar**2+eps_r**2) & + - eps_r*sin(theta)**2/(1._DP-eps_r**2) & + ) / ((1._DP + eps_r*cos(theta))**2) & + * sqrt(q_bar**2+eps_r**2) / q_bar / r_major + gdomgdy = 0._DP + ggxx = (q_0/q_bar)**2 + ggxy = ( s_hat*wzz*q_0/q_bar - eps_r*sin(wzz)/(1._DP-eps_r**2) )*q_0/q_bar + ggxz = - sin(wzz)/(1._DP-eps_r**2)/r_major*q_0/q_bar + ggyy = (s_hat*wzz*q_0/q_bar)**2 & + - 2._DP*q_0/q_bar*s_hat*wzz*eps_r*sin(wzz)/(1._DP-eps_r**2) & + + (q_bar**2+eps_r**2)/((1._DP+eps_r*cos(theta))**2)/(q_0**2) & + + (eps_r*sin(wzz))**2/(1._DP-eps_r**2)**2 + ggyz = ( -s_hat*wzz*q_0/q_bar*sin(wzz)/(1._DP-eps_r**2) & + + ((q_bar/q_0)**2)/((1._DP+eps_r*cos(theta))**2)/eps_r & + + eps_r*(sin(wzz)**2)/((1._DP-eps_r**2)**2) & + ) / r_major + ggzz = ( ((q_bar/q_0)**2)/((1._DP+eps_r*cos(theta))**2)/(eps_r**2) & + + (sin(wzz)**2)/((1._DP-eps_r**2)**2) & + ) / (r_major**2) + grootg_xyz = q_0*r_major*( 1._DP+eps_r*cos(theta) )**2 + + + else if( trim(equib_type) == "vmec" ) then + !- VMEC-BoozXform interface for stellarator equilibirum - + ! References on the previous implementation by VMEC-NEWBOZ is + ! [Ref.1] M. Nunami, T.-H. Watanabe, H. Sugama, Plasma Fusion Res. 5, + ! 016 (2010). + ! New interface for VMEC-BoozXform is developed by M. Nakata and + ! M. Nunami (Aug. 2016) in the same manner for IGS. + ! + ! Consider flux coordinates (rho,theta,zeta). + ! Using the toroidal flux psi_t, the normalized minor radius is + ! rho= sqrt(psi_t/psi_ta), and the minor radius at the last closed + ! flux surface is a=sqrt(2*psi_ta/B_ax). + ! Poloidal and toroidal angles are defined in the Boozer coordinates. + ! GKV coordinates (x,y,z) are + ! x = cx*(rho-rho0) + ! y = cy*(q(r)*theta-zeta) + ! z = theta + ! with cx=a, cy=cx*rho0/q0. + ! In these definitions, the factor on the magnetic field + ! B = cb * \nabla x \times \nabla y is + ! cb = (psi_p(rho))'/(cx*cy) = B_ax. + ! Normalized omg = B(z)/B_ax and cb = 1 in the B_ax unit. + ! The reference length is set to be r_major at the magnetic axis R0. + !- + r_major = 1._DP ! in the R0 unit + q_bar = q_0 + isw = 1 + call vmecbzx_boozx_coeff( isw, nss, ntheta, nzeta, s_input, iz, wzz, lz_l, & ! input + s_0, q_0, s_hat, eps_r, phi_ax, & ! output + gomg, grootg_xyz, gdomgdx, gdomgdz, gdomgdy, & + ggxx, ggxy, ggxz, ggyy, ggyz, ggzz ) + ! NOTE: phi_ax axisymmetric toroidal angle is stored for vmec, rather than theta + theta = phi_ax + cx = eps_r/s_0 ! = eps_a = a/R0 + cy = cx*s_0/q_0 + cb = 1._DP + + else if( trim(equib_type) == "eqdsk" ) then + !- EQDSK-IGS interface for tokamak equilibirum - + ! [Ref.] M. Nakata, A. Matsuyama, N. Aiba, S. Maeyama, M. Nunami, + ! and T.-H. Watanabe, Plasma Fusion Res. 9, 1403029 (2014). + ! + ! Consider flux coordinates (rho,theta,zeta). + ! GKV coordinates (x,y,z) are + ! x = cx*(rho-rho0) + ! y = cy*(q(r)*theta-zeta) + ! z = theta + ! with cx=a, cy=cx*rho0/q0. All explanation is the same as that in + ! equib_type == "vmec", except that poloidal and toroidal angles have + ! a choice of freedom: Hamada, Boozer, or axisymmetric coordinates. + !- + r_major = 1._DP ! in the R0 unit + q_bar = q_0 + isw = 1 + call igs_coeff( isw, mc_type, nss, ntheta, s_input, wzz, lz_l, & ! input + s_0, q_0, s_hat, eps_r, theta, & ! output + gomg, grootg_xyz, gdomgdx, gdomgdz, gdomgdy, & + ggxx, ggxy, ggxz, ggyy, ggyz, ggzz ) + cx = eps_r/s_0 ! = eps_a = a/R0 + cy = cx*s_0/q_0 + cb = 1._DP + +!sakano_ring-dipole st 202303 + else if( trim(equib_type) == "ring" ) then + !- Ring dipole geometry - + ! [Ref.] J. Sakano, Master thesis, Nagoya University (in Japanese). + ! + ! Consider flux coordinates (Psi,Theta,phi), where the magnetic + ! poloidal flux Psi<0, the geometrical poloidal angle Theta = arctan(Z/(R-a)), + ! the azimuthal angle of the cylindrical coordinate phi. + ! There is a ring current in direction of phi at R=a. The field line + ! passing through (R,Z)=(R0,0) is picked up as a flux-tube domain. + ! + ! GKV coordinates (x,y,z) are (right-handed system) + ! x = cx*(Psi0 - Psi)/Psi0 + ! y = cy*phi + ! z = Theta + ! with cx=R0, cy=R0. Note that Psi0 is the magnetic poloidal flux + ! at the center of the considered flux-tube domain. + ! In these definitions, the factor on the magnetic field + ! B = cb * \nabla x \times \nabla y is cb = Psi0/(R0*R0) = B0, + ! where B0 is the magnetic field strength at (R,Z)=(R0,0). + ! Normalized omg = B(z)/B0 and cb = 1 in the B0 unit. + ! The reference length is set to be R0 (not the ring current at R=a). + ! The normalized parameter to specify the flux-tube is + ! ring_a = a / R0 + !- + r_major = 1._DP ! in the R0 unit + q_bar = 0._DP + theta = wzz + call ring_coordinates( ring_a, wzz, & ! input + gomg, ub_dot_grdb, ub_crs_grdb, ggxx, ggxy, & ! output + ggxz, ggyy, ggyz, ggzz, grootg_xyz, gdomgdx, gdomgdz ) + gdomgdy = 0._DP + cx = 1._DP + cy = 1._DP + cb = 1._DP +!sakano_ring-dipole end 202303 + + else + + write( olog, * ) " # wrong choice of the equilibrium " + call flush(olog) + call MPI_Finalize(ierr_mpi) + stop + + end if + + call mtr_global%init(iz, wzz, theta, gomg, & + gdomgdx, gdomgdy, gdomgdz, ggxx, ggxy, & + ggxz, ggyy, ggyz, ggzz, grootg_xyz, & + gdomgdr, gdomgdt, gdomgdq, ggrr, ggrt, & + ggrq, ggtt, ggtq, ggqq, grootg_rtq) + + end do ! iz loop ends + + call mtr_global%xyz2rtq + call mtr_fourier%init + call mtr_fourier%dft_rtq2coef(mtr_global) + + call mtr_local%copy_global(mtr_global) + + if ( rankg == 0 ) then + do iz = -global_nz, global_nz-1 + write( omtr, fmt="(f15.8,SP,256E24.14e3)") & + mtr_global%zz(iz), mtr_global%theta(iz), & + mtr_global%omg(iz), mtr_global%domgdx(iz), & + mtr_global%domgdy(iz), mtr_global%domgdz(iz), & + mtr_global%gxx(iz), mtr_global%gxy(iz), & + mtr_global%gxz(iz), mtr_global%gyy(iz), & + mtr_global%gyz(iz), mtr_global%gzz(iz), & + mtr_global%rootg_xyz(iz) + end do + !call flush(omtr) + do iz = -global_nz, global_nz-1 + write( omtf, fmt="(f15.8,SP,256E24.14e3)") & + mtr_global%zz(iz), mtr_global%theta(iz), & + mtr_global%omg(iz), mtr_global%domgdr(iz), & + mtr_global%domgdt(iz), mtr_global%domgdq(iz), & + mtr_global%grr(iz), mtr_global%grt(iz), & + mtr_global%grq(iz), mtr_global%gtt(iz), & + mtr_global%gtq(iz), mtr_global%gqq(iz), & + mtr_global%rootg_rtq(iz) + end do + !call flush(omtr) + end if + + !%%% For debug %%% + ! do iz = -nz, nz-1 + ! write( 990000000+rankg, fmt="(f15.8,SP,256E24.14e3)") & + ! mtr_local%zz(iz), mtr_local%theta(iz), & + ! mtr_local%omg(iz), mtr_local%domgdx(iz), & + ! mtr_local%domgdy(iz), mtr_local%domgdz(iz), & + ! mtr_local%gxx(iz), mtr_local%gxy(iz), & + ! mtr_local%gxz(iz), mtr_local%gyy(iz), & + ! mtr_local%gyz(iz), mtr_local%gzz(iz), & + ! mtr_local%rootg_xyz(iz) + ! end do + ! call mtr_global%rtq2xyz + ! call mtr_global%xyz2rtq + ! call mtr_local%init(mtr_fourier, time_shearflow=0._DP) + ! if ( rankg == 0 ) then + ! do iz = -global_nz, global_nz-1 + ! write( 900000011, fmt="(f15.8,SP,256E24.14e3)") & + ! mtr_global%zz(iz), mtr_global%theta(iz), & + ! mtr_global%omg(iz), mtr_global%domgdx(iz), & + ! mtr_global%domgdy(iz), mtr_global%domgdz(iz), & + ! mtr_global%gxx(iz), mtr_global%gxy(iz), & + ! mtr_global%gxz(iz), mtr_global%gyy(iz), & + ! mtr_global%gyz(iz), mtr_global%gzz(iz), & + ! mtr_global%rootg_xyz(iz) + ! end do + ! do iz = -global_nz, global_nz-1 + ! write( 900000012, fmt="(f15.8,SP,256E24.14e3)") & + ! mtr_global%zz(iz), mtr_global%theta(iz), & + ! mtr_global%omg(iz), mtr_global%domgdr(iz), & + ! mtr_global%domgdt(iz), mtr_global%domgdq(iz), & + ! mtr_global%grr(iz), mtr_global%grt(iz), & + ! mtr_global%grq(iz), mtr_global%gtt(iz), & + ! mtr_global%gtq(iz), mtr_global%gqq(iz), & + ! mtr_global%rootg_rtq(iz) + ! end do + ! end if + ! do iz = -nz, nz-1 + ! write( 980000000+rankg, fmt="(f15.8,SP,256E24.14e3)") & + ! mtr_local%zz(iz), mtr_local%theta(iz), & + ! mtr_local%omg(iz), mtr_local%domgdx(iz), & + ! mtr_local%domgdy(iz), mtr_local%domgdz(iz), & + ! mtr_local%gxx(iz), mtr_local%gxy(iz), & + ! mtr_local%gxz(iz), mtr_local%gyy(iz), & + ! mtr_local%gyz(iz), mtr_local%gzz(iz), & + ! mtr_local%rootg_xyz(iz) + ! end do + ! do iz = -nz, nz-1 + ! write( 970000000+rankg, fmt="(f15.8,SP,256E24.14e3)") & + ! mtr_local%zz_labframe(iz), mtr_local%theta(iz), & + ! mtr_local%omg(iz), mtr_local%domgdr(iz), & + ! mtr_local%domgdt(iz), mtr_local%domgdq(iz), & + ! mtr_local%grr(iz), mtr_local%grt(iz), & + ! mtr_local%grq(iz), mtr_local%gtt(iz), & + ! mtr_local%gtq(iz), mtr_local%gqq(iz), & + ! mtr_local%rootg_rtq(iz) + ! end do + !%%%%%%%%%%%%%%%%%% + + END SUBROUTINE geom_init_metric + + +!-------------------------------------- + SUBROUTINE geom_set_operators +!-------------------------------------- + implicit none + real(kind=DP) :: wzz ! The rotating flux tube coordinate (= z'') + real(kind=DP) :: zz_lab ! The flux-coordinate theta in the lab frame (= z''+t*gamma_e/s_hat) + real(kind=DP) :: kkx, kky, domgdz, domgdx, domgdy + real(kind=DP) :: bb, kmo + real(kind=DP) :: gg0 + + real(kind=DP) :: cfsrf_l + real(kind=DP), dimension(1:3,1:3) :: gg + complex(kind=DP), dimension(:,:,:,:,:), allocatable :: wf + complex(kind=DP), dimension(:,:,:), allocatable :: nw + real(kind=DP), dimension(:,:,:), allocatable :: ww + + integer :: mx, my, iz, iv, im, is + + do iz = -nz, nz-1 + + wzz = zz(iz) + zz_lab = mtr_local%zz_labframe(iz) + omg(iz) = mtr_local%omg(iz) + domgdx = mtr_local%domgdx(iz) + domgdy = mtr_local%domgdy(iz) + domgdz = mtr_local%domgdz(iz) + gg(1,1) = mtr_local%gxx(iz) + gg(1,2) = mtr_local%gxy(iz) + gg(1,3) = mtr_local%gxz(iz) + gg(2,2) = mtr_local%gyy(iz) + gg(2,3) = mtr_local%gyz(iz) + gg(3,3) = mtr_local%gzz(iz) + rootg(iz) = mtr_local%rootg_xyz(iz) + gg(2,1) = gg(1,2) + gg(3,1) = gg(1,3) + gg(3,2) = gg(2,3) + +!!! for slab model + if ( trim(equib_type) == "slab") then + + dpara(iz) = dz * q_0 * r_major + + do im = 0, nm + vp(iz,im) = sqrt( 2._DP * mu(im) )!* omg(iz) ) + mir(iz,im) = 0._DP + do iv = 1, 2*nv + vdx(iz,iv,im) = 0._DP + vdy(iz,iv,im) = 0._DP + 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 ) ) + end do + end do ! im loop ends + + ksq(:,:,iz) = 0._DP + do my = ist_y, iend_y + do mx = -nx, nx + ksq(mx,my,iz) = kx(mx)**2 + ky(my)**2 + end do + end do + +!!! for the concentric and large-aspect-ratio model !!! + else if( trim(equib_type) == "analytic" ) then + + dpara(iz) = dz * q_0 * r_major + + 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 ) & + * ( sin(zz_lab) & + + eps_hor * lmmq * sin( lmmq * zz_lab - malpha ) & + + eps_mor * lmmqm1 * sin( lmmqm1 * zz_lab - malpha ) & + + eps_por * lmmqp1 * sin( lmmqp1 * zz_lab - malpha ) ) + + do iv = 1, 2*nv + vdx(iz,iv,im)= & + - ( vl(iv)**2 + omg(iz)*mu(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 ) & + + rdeps3_10 * cos( lmmqp1 * zz_lab - malpha ) ) & + + ( 1._DP + s_hat * wzz * 0._DP ) & + * ( sin( zz_lab ) & + + eps_hor * lprd * sin( lmmq * zz_lab - malpha ) & + + eps_mor * lprdm1 * sin( lmmqm1 * zz_lab - malpha ) & + + eps_por * lprdp1 * sin( lmmqp1 * zz_lab - malpha ) ) & + ) * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) + + vdy(iz,iv,im)= & + - ( vl(iv)**2 + omg(iz)*mu(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 ) & + + rdeps3_10 * cos( lmmqp1 * zz_lab - malpha ) ) & + + ( 0._DP + s_hat * wzz * 1._DP ) & + * ( sin( zz_lab ) & + + eps_hor * lprd * sin( lmmq * zz_lab - malpha ) & + + eps_mor * lprdm1 * sin( lmmqm1 * zz_lab - malpha ) & + + eps_por * lprdp1 * sin( lmmqp1 * zz_lab - malpha ) ) & + ) * ( 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 ) ) + + end do + + end do ! im loop ends + + ksq(:,:,iz) = 0._DP + do my = ist_y, iend_y + do mx = -nx, nx + ksq(mx,my,iz) = ( kx(mx) + s_hat * wzz * ky(my) )**2 + ky(my)**2 + end do + end do + +!!! for s-alpha !!! <--- the current version is the same as "analytic" + else if( trim(equib_type) == "s-alpha" .or. trim(equib_type) == "s-alpha-shift" ) then + + dpara(iz) = dz* q_0 * r_major + + kkx = -r_major * (q_0/q_bar) & + * ( gg(1,1)*gg(2,3) - gg(1,2)*gg(1,3) )*domgdz + kky = r_major * (q_bar/q_0) & + * ( domgdx - ( gg(1,2)*gg(2,3) - gg(2,2)*gg(1,3) )*domgdz ) + + do im = 0, nm + + vp(iz,im) = sqrt( 2._DP * mu(im) * omg(iz) ) + + mir(iz,im) = mu(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 & + * kkx & + * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) + + vdy(iz,iv,im) = & + ( vl(iv)**2 + omg(iz)*mu(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 ) ) & + * (q_bar/q_0) + end do + + end do ! im loop ends + + + ksq(:,:,iz) = 0._DP + do my = ist_y, iend_y + do mx = -nx, nx + ksq(mx,my,iz) = ( kx(mx) + ( s_hat * wzz - alpha_MHD*sin(zz_lab) ) & + * ky(my) )**2 + ky(my)**2 ! with Shafranov shift + end do + end do + +!!! for circular MHD equilibrium !!! + else if( trim(equib_type) == "circ-MHD" ) then + + dpara(iz) = dz * omg(iz) * rootg(iz) + + kkx = r_major*( -domgdy + (gg(1,3)*gg(1,2) - gg(1,1)*gg(2,3))*domgdz/omg(iz)**2 ) + kky = r_major*( domgdx + (gg(1,3)*gg(2,2) - gg(1,2)*gg(2,3))*domgdz/omg(iz)**2 ) + + do im = 0, nm + + vp(iz,im) = sqrt( 2._DP * mu(im) * omg(iz) ) + + mir(iz,im) = mu(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) ) & + * kkx & + * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) + + vdy(iz,iv,im)= & + ( vl(iv)**2 + omg(iz)*mu(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 ) ) + end do + + end do ! im loop ends + + ksq(:,:,iz) = 0._DP + do my = ist_y, iend_y + do mx = -nx, nx + ksq(mx,my,iz) = (kx(mx)**2)*gg(1,1) & + + 2._DP*kx(mx)*ky(my)*gg(1,2) & + + (ky(my)**2)*gg(2,2) + end do + end do + +! this is new vmec-BoozXform interface by M. Nakata & M. Nunami (Aug. 2016) + else if( trim(equib_type) == "vmec" ) then + + dpara(iz) = dz * omg(iz) * rootg(iz) + + kkx = r_major*( -domgdy + (gg(1,3)*gg(1,2) - gg(1,1)*gg(2,3))*domgdz/omg(iz)**2 ) + kky = r_major*( domgdx + (gg(1,3)*gg(2,2) - gg(1,2)*gg(2,3))*domgdz/omg(iz)**2 ) + + do im = 0, nm + + vp(iz,im) = sqrt( 2._DP * mu(im) * omg(iz) ) + + mir(iz,im) = mu(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) ) & + * kkx & + * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) + + + vdy(iz,iv,im) = & + ( vl(iv)**2 + omg(iz)*mu(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 + * ( beta*(R0_Ln(ranks) + R0_Lt(ranks)) ) & + * ( 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 ) ) + end do + + end do ! im loop ends + + ksq(:,:,iz) = 0._DP + do my = ist_y, iend_y + do mx = -nx, nx + ksq(mx,my,iz) = (kx(mx)**2)*gg(1,1) & + + 2._DP*kx(mx)*ky(my)*gg(1,2) & + + (ky(my)**2)*gg(2,2) + end do + end do + + else if( trim(equib_type) == "eqdsk" ) then + + dpara(iz) = dz * omg(iz) * rootg(iz) + + kkx = r_major*( -domgdy + (gg(1,3)*gg(1,2) - gg(1,1)*gg(2,3))*domgdz/omg(iz)**2 ) + kky = r_major*( domgdx + (gg(1,3)*gg(2,2) - gg(1,2)*gg(2,3))*domgdz/omg(iz)**2 ) + + do im = 0, nm + + vp(iz,im) = sqrt( 2._DP * mu(im) * omg(iz) ) + + mir(iz,im) = mu(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) ) & + * kkx & + * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) + + vdy(iz,iv,im) = & + ( vl(iv)**2 + omg(iz)*mu(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 + * ( beta*(R0_Ln(ranks) + R0_Lt(ranks)) ) & + * ( 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 ) ) + end do + + end do ! im loop ends + + ksq(:,:,iz) = 0._DP + do my = ist_y, iend_y + do mx = -nx, nx + ksq(mx,my,iz) = (kx(mx)**2)*gg(1,1) & + + 2._DP*kx(mx)*ky(my)*gg(1,2) & + + (ky(my)**2)*gg(2,2) + end do + end do + +!sakano_ring-dipole st 202303 + else if( trim(equib_type) == "ring" ) then + + dpara(iz) = dz * omg(iz) * rootg(iz) + + kkx = 0._DP + kky = r_major*( domgdx + (gg(1,3)*gg(2,2) - gg(1,2)*gg(2,3))*domgdz/omg(iz)**2 ) + + do im = 0, nm +! 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) ) + + !mir(iz,im) = mu(im) * ub_dot_grdb + mir(iz,im) = mu(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) ) & + ! * ( 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) ) & + * 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 + end do + + end do ! im loop ends + + ksq(:,:,iz) = 0._DP + do my = ist_y, iend_y + do mx = -nx, nx + ksq(mx,my,iz) = ( kx(mx) * gg(1,1) )**2 + ( ky(my) * gg(2,2) )**2 + end do + end do +!sakano_ring-dipole end 202303 + + else + + write( olog, * ) " # wrong choice of the equilibrium " + call flush(olog) + call MPI_Finalize(ierr_mpi) + stop + + end if + + + 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) ) & + * 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) ) + call math_j2( kmo, j2(mx,my,iz,im) ) + end do + end do + end do + + + do my = ist_y, iend_y + do mx = -nx, nx + bb = ksq(mx,my,iz) / omg(iz)**2 & + * tau(ranks)*Anum(ranks)/(Znum(ranks)**2) + call math_g0( bb, g0(mx,my,iz) ) + end do + end do + + end do ! iz loop ends + + cfsrf = 0._DP + cfsrf_l = 0._DP + do iz = -nz, nz-1 + cfsrf_l = cfsrf_l + rootg(iz) + ! normalization coefficient for + ! the surface average + end do + call MPI_Allreduce( cfsrf_l, cfsrf, 1, MPI_DOUBLE_PRECISION, & + MPI_SUM, zsp_comm_world, ierr_mpi ) + + if ( vel_rank == 0 ) then + do iz = -nz, nz-1 + dvp(iz) = sqrt( 2._DP * (0.5_DP * dm**2) * omg(iz) ) + end do + end if + call MPI_Bcast( dvp, 2*nz, MPI_DOUBLE_PRECISION, 0, & + vel_comm_world, ierr_mpi ) + + 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) ) & + / sqrt( twopi**3 ) + end do + end do + end do + + allocate( ww(-nx:nx,0:ny,-nz:nz-1) ) + +! --- GK polarization factor for efield calculation + fct_poisson(:,:,:) = 0._DP + fct_e_energy(:,:,:) = 0._DP + + ww(:,:,:) = 0._DP + do iz = -nz, nz-1 + do my = ist_y, iend_y + do mx = -nx, nx + + if ( rankw == 0 .and. mx == 0 .and. my == 0 ) then !- (0,0) mode + + fct_poisson(mx,my,iz) = 0._DP + fct_e_energy(mx,my,iz) = 0._DP + + else + + ww(mx,my,iz) = lambda_i * ksq(mx,my,iz) + do is = 0, ns-1 + bb = ksq(mx,my,iz) / omg(iz)**2 & + * tau(is)*Anum(is)/(Znum(is)**2) + call math_g0( bb, gg0 ) + ww(mx,my,iz) = ww(mx,my,iz) & + + Znum(is) * fcs(is) / tau(is) * ( 1._DP - gg0 ) + end do + fct_poisson(mx,my,iz) = 1._DP / ww(mx,my,iz) + fct_e_energy(mx,my,iz) = ww(mx,my,iz) + + end if + + end do + end do + end do + + +! --- ZF-factor for adiabatic model + if ( ns == 1 ) then + + ww(:,:,:) = 0._DP + do iz = -nz, nz-1 + my = 0 + do mx = -nx, nx + ww(mx,my,iz) = ( 1._DP - g0(mx,my,iz) ) & + / ( 1._DP - g0(mx,my,iz) + tau(0)*tau_ad ) + end do + end do + + call intgrl_fsrf ( ww, fctgt ) + + if ( rankw == 0 ) then + fctgt(0) = ( 1._DP - g0(0,0,0) ) / ( 1._DP - g0(0,0,0) + tau(0)*tau_ad ) + ! g0(0,0,iz) has no z dependence + endif + + endif + + deallocate( ww ) + + allocate( wf(-nx:nx,0:ny,-nz:nz-1,1:2*nv,0:nm) ) + allocate( nw(-nx:nx,0:ny,-nz:nz-1) ) + wf(:,:,:,:,:) = ( 0._DP, 0._DP ) + nw(:,:,:) = ( 0._DP, 0._DP ) + +! --- GK polarization factor for mfield calculation + fct_ampere(:,:,:) = 0._DP + fct_m_energy(:,:,:) = 0._DP + + if ( beta .ne. 0._DP ) 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 + wf(mx,my,iz,iv,im) = Znum(ranks) * fcs(ranks) / Anum(ranks) & + * vl(iv)**2 * j0(mx,my,iz,im)**2 * fmx(iz,iv,im) + end do + end do + end do + end do + end do + + call intgrl_v0_moment_ms ( wf, nw ) + + do iz = -nz, nz-1 + do my = ist_y, iend_y + do mx = -nx, nx + fct_ampere(mx,my,iz) = 1._DP / real( ksq(mx,my,iz) + beta * nw(mx,my,iz), kind=DP ) + fct_m_energy(mx,my,iz) = ksq(mx,my,iz) / beta + end do + end do + end do + + if ( rankw == 0 ) then + do iz = -nz, nz-1 + fct_ampere(0,0,iz) = 0._DP + fct_m_energy(0,0,iz) = 0._DP + end do + end if + + end if + + deallocate( wf ) + deallocate( nw ) + + END SUBROUTINE geom_set_operators + +!-------------------------------------- + SUBROUTINE geom_reset_time(time_shearflow) +!-------------------------------------- + implicit none + real(kind=DP), intent(in) :: time_shearflow + call mtr_local%init(mtr_fourier, time_shearflow) + call geom_set_operators + !NOTE: colliimp_set_param in GKV_colliimp should also be updated. + END SUBROUTINE geom_reset_time + +!-------------------------------------- + SUBROUTINE geom_increment_time(dt_shearflow) +!-------------------------------------- + implicit none + real(kind=DP), intent(in) :: dt_shearflow + call mtr_local%update(mtr_fourier, dt_shearflow) + call geom_set_operators + !NOTE: colliimp_set_param in GKV_colliimp should also be updated. + END SUBROUTINE geom_increment_time + + +!-------------------------------------- + SUBROUTINE metric_global_init(self, iz, wzz, theta, gomg, & + gdomgdx, gdomgdy, gdomgdz, ggxx, ggxy, & + ggxz, ggyy, ggyz, ggzz, grootg_xyz, & + gdomgdr, gdomgdt, gdomgdq, ggrr, ggrt, & + ggrq, ggtt, ggtq, ggqq, grootg_rtq) +!-------------------------------------- + implicit none + class(metric_global), intent(inout) :: self + integer, intent(in) :: iz + real(kind=DP), intent(in) :: wzz, theta, gomg + real(kind=DP), intent(in) :: gdomgdx, gdomgdy, gdomgdz, & + ggxx, ggxy, ggxz, ggyy, ggyz, ggzz, grootg_xyz + real(kind=DP), intent(in) :: gdomgdr, gdomgdt, gdomgdq, & + ggrr, ggrt, ggrq, ggtt, ggtq, ggqq, grootg_rtq + + self%zz(iz) = wzz + self%theta(iz) = theta + self%omg(iz) = gomg + self%domgdx(iz) = gdomgdx + self%domgdy(iz) = gdomgdy + self%domgdz(iz) = gdomgdz + self%gxx(iz) = ggxx + self%gxy(iz) = ggxy + self%gxz(iz) = ggxz + self%gyy(iz) = ggyy + self%gyz(iz) = ggyz + self%gzz(iz) = ggzz + self%rootg_xyz(iz) = grootg_xyz + self%domgdr(iz) = gdomgdr + self%domgdt(iz) = gdomgdt + self%domgdq(iz) = gdomgdq + self%grr(iz) = ggrr + self%grt(iz) = ggrt + self%grq(iz) = ggrq + self%gtt(iz) = ggtt + self%gtq(iz) = ggtq + self%gqq(iz) = ggqq + self%rootg_rtq(iz) = grootg_rtq + + END SUBROUTINE metric_global_init + + +!-------------------------------------- + SUBROUTINE metric_global_xyz2rtq(self) +!-------------------------------------- + implicit none + class(metric_global), intent(inout) :: self + real(kind=DP) :: wzz + real(kind=DP) :: gdomgdx, gdomgdy, gdomgdz, & + ggxx, ggxy, ggxz, ggyy, ggyz, ggzz, grootg_xyz + real(kind=DP) :: gdomgdr, gdomgdt, gdomgdq, & + ggrr, ggrt, ggrq, ggtt, ggtq, ggqq, grootg_rtq + integer :: iz + + do iz = -global_nz, global_nz-1 + ! load (x,y,z) + wzz = self%zz(iz) + gdomgdx = self%domgdx(iz) + gdomgdy = self%domgdy(iz) + gdomgdz = self%domgdz(iz) + ggxx = self%gxx(iz) + ggxy = self%gxy(iz) + ggxz = self%gxz(iz) + ggyy = self%gyy(iz) + ggyz = self%gyz(iz) + ggzz = self%gzz(iz) + grootg_xyz = self%rootg_xyz(iz) + + ! translate (x,y,z)->(r,t,q)=(rho,theta,zeta) + ! NOTE: cx*rho0/(cy*q_0=1) is used. + gdomgdr = cx*gdomgdx + cx*s_hat*wzz*gdomgdy + gdomgdt = gdomgdz + cy*q_0*gdomgdy + gdomgdq = - cy*gdomgdy + ggrr = ggxx/cx**2 + ggrt = ggxz/cx + ggrq = (s_hat*wzz*ggxx-ggxy)/(cx*cy) + q_0*ggxz/cx + ggtt = ggzz + ggtq = (s_hat*wzz*ggxz-ggyz)/cy + q_0*ggzz + ggqq = (s_hat*wzz/cy)**2*ggxx - 2._DP*(s_hat*wzz/cy**2)*ggxy & + + 2._DP*(q_0*s_hat*wzz/cy)*ggxz + ggyy/cy**2 & + - 2._DP*(q_0/cy)*ggyz + q_0**2*ggzz + grootg_rtq = cx*cy*grootg_xyz + + ! store (r,t,q) + self%domgdr(iz) = gdomgdr + self%domgdt(iz) = gdomgdt + self%domgdq(iz) = gdomgdq + self%grr(iz) = ggrr + self%grt(iz) = ggrt + self%grq(iz) = ggrq + self%gtt(iz) = ggtt + self%gtq(iz) = ggtq + self%gqq(iz) = ggqq + self%rootg_rtq(iz) = grootg_rtq + end do + + END SUBROUTINE metric_global_xyz2rtq + + +!-------------------------------------- + SUBROUTINE metric_global_rtq2xyz(self) +!-------------------------------------- + implicit none + class(metric_global), intent(inout) :: self + real(kind=DP) :: wzz + real(kind=DP) :: gdomgdx, gdomgdy, gdomgdz, & + ggxx, ggxy, ggxz, ggyy, ggyz, ggzz, grootg_xyz + real(kind=DP) :: gdomgdr, gdomgdt, gdomgdq, & + ggrr, ggrt, ggrq, ggtt, ggtq, ggqq, grootg_rtq + integer :: iz + + do iz = -global_nz, global_nz-1 + ! load (r,t,q)=(rho,theta,zeta) + wzz = self%zz(iz) + gdomgdr = self%domgdr(iz) + gdomgdt = self%domgdt(iz) + gdomgdq = self%domgdq(iz) + ggrr = self%grr(iz) + ggrt = self%grt(iz) + ggrq = self%grq(iz) + ggtt = self%gtt(iz) + ggtq = self%gtq(iz) + ggqq = self%gqq(iz) + grootg_rtq = self%rootg_rtq(iz) + + ! translate (r,t,q)->(x,y,z) + ! NOTE: cx*rho0/(cy*q_0=1) is used. + gdomgdx = gdomgdr/cx + s_hat*wzz*gdomgdq/cy + gdomgdy = - gdomgdq/cy + gdomgdz = gdomgdt + q_0*gdomgdq + ggxx = cx**2*ggrr + ggxy = cx**2*s_hat*wzz*ggrr + cx*cy*(q_0*ggrt - ggrq) + ggxz = cx*ggrt + ggyy = (cx*s_hat*wzz)**2*ggrr + 2._DP*cx*cy*s_hat*wzz*(q_0*ggrt-ggrq) & + + (cy*q_0)**2*ggtt - 2._DP*cy**2*q_0*ggtq + cy**2*ggqq + ggyz = cx*s_hat*wzz*ggrt + cy*q_0*ggtt - cy*ggtq + ggzz = ggtt + grootg_xyz = grootg_rtq/(cx*cy) + + ! store (x,y,z) + self%domgdx(iz) = gdomgdx + self%domgdy(iz) = gdomgdy + self%domgdz(iz) = gdomgdz + self%gxx(iz) = ggxx + self%gxy(iz) = ggxy + self%gxz(iz) = ggxz + self%gyy(iz) = ggyy + self%gyz(iz) = ggyz + self%gzz(iz) = ggzz + self%rootg_xyz(iz) = grootg_xyz + end do + + END SUBROUTINE metric_global_rtq2xyz + + +!-------------------------------------- + SUBROUTINE metric_fourier_init(self) +!-------------------------------------- + implicit none + class(metric_fourier), intent(inout) :: self + real(kind=DP) :: kzmin + integer :: iz + + kzmin = 2._DP * pi / (2._DP * lz) + do iz = -global_nz, global_nz-1 + self%kz(iz) = iz * kzmin + end do + + END SUBROUTINE metric_fourier_init + + +!-------------------------------------- + SUBROUTINE forward_dft_globalz(zz_global,kz,fz,fk) +!-------------------------------------- + implicit none + real(kind=DP), intent(in), & + dimension(-global_nz:global_nz-1) :: zz_global, kz, fz + complex(kind=DP), intent(out), & + dimension(-global_nz:global_nz-1) :: fk + integer :: iz, mz + + fk(:) = (0._DP, 0._DP) + do mz = -global_nz, global_nz-1 + do iz = -global_nz, global_nz-1 + fk(mz) = fk(mz) + fz(iz)*exp(-ui*kz(mz)*zz_global(iz))*dz/(2._DP*lz) + end do + end do + + END SUBROUTINE forward_dft_globalz + + +!!-------------------------------------- +! SUBROUTINE backward_dft_globalz(zz_global,kz,fk,fz) +!!-------------------------------------- +! implicit none +! real(kind=DP), intent(in), & +! dimension(-global_nz:global_nz-1) :: zz_global, kz +! complex(kind=DP), intent(in), & +! dimension(-global_nz:global_nz-1) :: fk +! real(kind=DP), intent(out), & +! dimension(-global_nz:global_nz-1) :: fz +! integer :: iz, mz +! +! fz(:) = 0._DP +! do iz = -global_nz, global_nz-1 +! do mz = -global_nz, global_nz-1 +! fz(iz) = fz(iz) + real(fk(mz)*exp(ui*kz(mz)*zz_global(iz)), kind=DP) +! end do +! end do +! +! END SUBROUTINE backward_dft_globalz + + +!-------------------------------------- + SUBROUTINE backward_dft_localz(zz_local,kz,fk,fz) +!-------------------------------------- + implicit none + real(kind=DP), intent(in), & + dimension(-nz:nz-1) :: zz_local + real(kind=DP), intent(in), & + dimension(-global_nz:global_nz-1) :: kz + complex(kind=DP), intent(in), & + dimension(-global_nz:global_nz-1) :: fk + real(kind=DP), intent(out), & + dimension(-nz:nz-1) :: fz + integer :: iz, mz + + fz(:) = 0._DP + do iz = -nz, nz-1 + do mz = -global_nz, global_nz-1 + fz(iz) = fz(iz) + real(fk(mz)*exp(ui*kz(mz)*zz_local(iz)), kind=DP) + end do + end do + + END SUBROUTINE backward_dft_localz + + +!-------------------------------------- + SUBROUTINE metric_fourier_dft_rtq2coef(self, mtr_g) +!-------------------------------------- + implicit none + class(metric_fourier), intent(inout) :: self + class(metric_global), intent(in) :: mtr_g + real(kind=DP), dimension(-global_nz:global_nz-1) :: theta_tilde + + ! theta = zz + theta_tilde(zz), theta_tilde is a periodic function. + if (trim(equib_type) == "vmec") then + theta_tilde = mtr_g%theta - q_0 * mtr_g%zz ! Axisymmetric toroidal angle phi_ax + else + theta_tilde = mtr_g%theta - mtr_g%zz + end if + call forward_dft_globalz(mtr_g%zz, self%kz, theta_tilde, self%theta_tilde) + call forward_dft_globalz(mtr_g%zz, self%kz, mtr_g%omg , self%omg ) + call forward_dft_globalz(mtr_g%zz, self%kz, mtr_g%domgdr , self%domgdr ) + call forward_dft_globalz(mtr_g%zz, self%kz, mtr_g%domgdt , self%domgdt ) + call forward_dft_globalz(mtr_g%zz, self%kz, mtr_g%domgdq , self%domgdq ) + call forward_dft_globalz(mtr_g%zz, self%kz, mtr_g%grr , self%grr ) + call forward_dft_globalz(mtr_g%zz, self%kz, mtr_g%grt , self%grt ) + call forward_dft_globalz(mtr_g%zz, self%kz, mtr_g%grq , self%grq ) + call forward_dft_globalz(mtr_g%zz, self%kz, mtr_g%gtt , self%gtt ) + call forward_dft_globalz(mtr_g%zz, self%kz, mtr_g%gtq , self%gtq ) + call forward_dft_globalz(mtr_g%zz, self%kz, mtr_g%gqq , self%gqq ) + call forward_dft_globalz(mtr_g%zz, self%kz, mtr_g%rootg_rtq, self%rootg_rtq) + ! NOTE: + ! Arguments are (zz_global(in),kz_global(in),omg_global(in),coef_global(out)) + + END SUBROUTINE metric_fourier_dft_rtq2coef + + +!-------------------------------------- + SUBROUTINE metric_local_dft_coef2rtq(self, mtr_f) +!-------------------------------------- + implicit none + class(metric_local), intent(inout) :: self + class(metric_fourier), intent(in) :: mtr_f + real(kind=DP), dimension(-nz:nz-1) :: theta_tilde + + ! theta = zz + theta_tilde(zz), theta_tilde is a periodic function. + call backward_dft_localz(self%zz_labframe, mtr_f%kz, mtr_f%theta_tilde, theta_tilde ) + if (trim(equib_type) == "vmec") then + self%theta = q_0 * self%zz_labframe + theta_tilde ! Axisymmetric toroidal angle phi_ax = q_0*zz + phi_tilde(zz) + else + self%theta = self%zz_labframe + theta_tilde + end if + call backward_dft_localz(self%zz_labframe, mtr_f%kz, mtr_f%omg , self%omg ) + call backward_dft_localz(self%zz_labframe, mtr_f%kz, mtr_f%domgdr , self%domgdr ) + call backward_dft_localz(self%zz_labframe, mtr_f%kz, mtr_f%domgdt , self%domgdt ) + call backward_dft_localz(self%zz_labframe, mtr_f%kz, mtr_f%domgdq , self%domgdq ) + call backward_dft_localz(self%zz_labframe, mtr_f%kz, mtr_f%grr , self%grr ) + call backward_dft_localz(self%zz_labframe, mtr_f%kz, mtr_f%grt , self%grt ) + call backward_dft_localz(self%zz_labframe, mtr_f%kz, mtr_f%grq , self%grq ) + call backward_dft_localz(self%zz_labframe, mtr_f%kz, mtr_f%gtt , self%gtt ) + call backward_dft_localz(self%zz_labframe, mtr_f%kz, mtr_f%gtq , self%gtq ) + call backward_dft_localz(self%zz_labframe, mtr_f%kz, mtr_f%gqq , self%gqq ) + call backward_dft_localz(self%zz_labframe, mtr_f%kz, mtr_f%rootg_rtq, self%rootg_rtq) + ! NOTE: + ! Arguments are (zz_local(in),kz_global(in), coef_global(in), omg_local(out)). + ! Fourier coefficients have been evaluated in the lab frame at t=0. + ! self%zz_labframe (= z''+t*gamma_e/s_hat) is the time-dependent flux-coordinate theta in the lab frame. + + END SUBROUTINE metric_local_dft_coef2rtq + + +!-------------------------------------- + SUBROUTINE metric_local_rtq2xyz(self) +!-------------------------------------- + implicit none + class(metric_local), intent(inout) :: self + real(kind=DP) :: wzz + real(kind=DP) :: gdomgdx, gdomgdy, gdomgdz, & + ggxx, ggxy, ggxz, ggyy, ggyz, ggzz, grootg_xyz + real(kind=DP) :: gdomgdr, gdomgdt, gdomgdq, & + ggrr, ggrt, ggrq, ggtt, ggtq, ggqq, grootg_rtq + integer :: iz + + do iz = -nz, nz-1 + ! load (r,t,q)=(rho,theta,zeta) + wzz = self%zz(iz) + gdomgdr = self%domgdr(iz) + gdomgdt = self%domgdt(iz) + gdomgdq = self%domgdq(iz) + ggrr = self%grr(iz) + ggrt = self%grt(iz) + ggrq = self%grq(iz) + ggtt = self%gtt(iz) + ggtq = self%gtq(iz) + ggqq = self%gqq(iz) + grootg_rtq = self%rootg_rtq(iz) + + ! translate (r,t,q)->(x,y,z) + ! NOTE: cx*rho0/(cy*q_0=1) is used. + gdomgdx = gdomgdr/cx + s_hat*wzz*gdomgdq/cy + gdomgdy = - gdomgdq/cy + gdomgdz = gdomgdt + q_0*gdomgdq + ggxx = cx**2*ggrr + ggxy = cx**2*s_hat*wzz*ggrr + cx*cy*(q_0*ggrt - ggrq) + ggxz = cx*ggrt + ggyy = (cx*s_hat*wzz)**2*ggrr + 2._DP*cx*cy*s_hat*wzz*(q_0*ggrt-ggrq) & + + (cy*q_0)**2*ggtt - 2._DP*cy**2*q_0*ggtq + cy**2*ggqq + ggyz = cx*s_hat*wzz*ggrt + cy*q_0*ggtt - cy*ggtq + ggzz = ggtt + grootg_xyz = grootg_rtq/(cx*cy) + + ! store (x,y,z) + self%domgdx(iz) = gdomgdx + self%domgdy(iz) = gdomgdy + self%domgdz(iz) = gdomgdz + self%gxx(iz) = ggxx + self%gxy(iz) = ggxy + self%gxz(iz) = ggxz + self%gyy(iz) = ggyy + self%gyz(iz) = ggyz + self%gzz(iz) = ggzz + self%rootg_xyz(iz) = grootg_xyz + end do + + END SUBROUTINE metric_local_rtq2xyz + + +!-------------------------------------- + SUBROUTINE metric_local_copy_global(self, mtr_g) +!-------------------------------------- + implicit none + class(metric_local), intent(inout) :: self + class(metric_global), intent(in) :: mtr_g + integer :: iz, giz + + do iz = -nz, nz-1 + giz = iz - global_nz + 2*nz * rankz + nz + self%zz_labframe(iz) = mtr_g%zz(giz) + self%zz(iz) = mtr_g%zz(giz) + self%theta(iz) = mtr_g%theta(giz) + self%omg(iz) = mtr_g%omg(giz) + self%domgdx(iz) = mtr_g%domgdx(giz) + self%domgdy(iz) = mtr_g%domgdy(giz) + self%domgdz(iz) = mtr_g%domgdz(giz) + self%gxx(iz) = mtr_g%gxx(giz) + self%gxy(iz) = mtr_g%gxy(giz) + self%gxz(iz) = mtr_g%gxz(giz) + self%gyy(iz) = mtr_g%gyy(giz) + self%gyz(iz) = mtr_g%gyz(giz) + self%gzz(iz) = mtr_g%gzz(giz) + self%rootg_xyz(iz) = mtr_g%rootg_xyz(giz) + end do + + END SUBROUTINE metric_local_copy_global + + +!-------------------------------------- + SUBROUTINE metric_local_init(self, mtr_f, time_shearflow) +!-------------------------------------- + implicit none + class(metric_local), intent(inout) :: self + class(metric_fourier), intent(in) :: mtr_f + real(kind=DP), intent(in) :: time_shearflow + + self%zz(:) = zz(:) + self%zz_labframe(:) = zz(:) + time_shearflow * gamma_e / s_hat + call self%dft_coef2rtq(mtr_f) + call self%rtq2xyz + + END SUBROUTINE metric_local_init + + +!-------------------------------------- + SUBROUTINE metric_local_update(self, mtr_f, dt_shearflow) +!-------------------------------------- + implicit none + class(metric_local), intent(inout) :: self + class(metric_fourier), intent(in) :: mtr_f + real(kind=DP), intent(in) :: dt_shearflow + + self%zz_labframe(:) = self%zz_labframe(:) + dt_shearflow * gamma_e / s_hat + call self%dft_coef2rtq(mtr_f) + call self%rtq2xyz + + END SUBROUTINE metric_local_update + + +END MODULE GKV_geom diff --git a/src/gkvp_header.f90 b/src/gkvp_header.f90 index ecf41e9..f40f1db 100644 --- a/src/gkvp_header.f90 +++ b/src/gkvp_header.f90 @@ -18,6 +18,11 @@ MODULE GKV_header ! ! Update history ! -------------- +! gkvp_f0.62 (S. Maeyama, Mar 2023) +! - flag_shearflow = "rotating" is set as a default. Alternatively, +! flag_shaerflow = "remap" is still available for time-discontinuous +! wave-vector remap with nearest grid approximation. +! - File I/O unit "omtf" is added for metrics in flux-coordinates. ! gkvp_f0.61 (S. Maeyama, Mar 2021) ! - equib_type is extended from len=8 to len=15. ! gkvp_f0.57 (S. Maeyama, Oct 2020) @@ -38,9 +43,9 @@ MODULE GKV_header ! in x, y,z,v,m (0:2*nxw-1, 0:2*nyw-1,-global_nz:global_nz-1,1:2*global_nv,0:global_nm) ! in kx,ky,z,v,m ( -nx:nx,0:global_ny,-global_nz:global_nz-1,1:2*global_nv,0:global_nm) - integer, parameter :: nxw = 2, nyw = 20 - integer, parameter :: nx = 0, global_ny = 12 ! 2/3 de-aliasing rule - integer, parameter :: global_nz = 48, global_nv = 24, global_nm = 15 + 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 :: nzb = 2, & ! the number of ghost grids in z nvb = 2 ! the number of ghost grids in v and m @@ -49,7 +54,7 @@ MODULE GKV_header ! Data distribution for MPI !-------------------------------------- - integer, parameter :: nprocw = 2, nprocz = 4, nprocv = 2, nprocm = 2, nprocs = 1 + integer, parameter :: nprocw = 1, nprocz = 2, nprocv = 4, nprocm = 2, nprocs = 1 !-------------------------------------- ! Parameters for variable sizes @@ -169,7 +174,7 @@ MODULE GKV_header tau, & ! T-ratio dns1 ! initial perturbation amp. real(kind=DP) :: dv, cfsrf, lambda_i, q_0, q_bar, beta, tau_ad, vmax - real(kind=DP) :: mach, uprime, gamma_e, kxmin_g, kymin_g, tlim_exb + real(kind=DP) :: mach, uprime, gamma_e, kxmin_g, kymin_g, tlim_exb, s_hat_g real(kind=DP) :: Nref, Lref, Tref, Zeff integer :: iFLR, icheck, ibprime, nx0 real(kind=DP) :: baxfactor @@ -201,6 +206,11 @@ MODULE GKV_header character(15) :: equib_type ! "analytic", "s-alpha", "s-alpha-shift", ! "circ-MHD", "vmec", "eqdsk", "slab" + !character(15) :: flag_shearflow = "remap" ! Wavevector remap method + ! ! with nearest grid approximation + ! ! (Discontinuous in time) + character(15) :: flag_shearflow = "rotating" ! Rotating flux tube model + ! --- unit numbers for I/O integer, parameter :: inml = 5, & olog = 10, & @@ -228,6 +238,7 @@ MODULE GKV_header inbz = 14, & ivmc = 15, & omtr = 16, & + omtf = 17, & ovmc = olog diff --git a/src/gkvp_main.f90 b/src/gkvp_main.f90 index 04e5db3..04c42e1 100644 --- a/src/gkvp_main.f90 +++ b/src/gkvp_main.f90 @@ -13,9 +13,9 @@ PROGRAM GKV_main ! | ! colli, colliimp, exb, shearflow ! | -! bndry, fft, fld, zfilter +! bndry, fft, fld, zfilter, geom ! | -! clock, intgrl, tips, freq, igs, vmecbzx, fileio +! clock, intgrl, tips, freq, igs, vmecbzx, ring, fileio ! | ! mpienv, math ! | @@ -136,7 +136,7 @@ PROGRAM GKV_main end if - if (gamma_e /= 0._DP) then + if (gamma_e /= 0._DP .and. trim(flag_shearflow) == "remap") then call shearflow_kxmap( time, ff, phi, Al, hh ) if (time > tlim_exb - eps .AND. cflg == 0 ) then write( olog, * ) "" diff --git a/src/gkvp_ring.f90 b/src/gkvp_ring.f90 new file mode 100644 index 0000000..7df18a1 --- /dev/null +++ b/src/gkvp_ring.f90 @@ -0,0 +1,474 @@ +MODULE ring_header + + implicit none + + integer, parameter :: DP = selected_real_kind(14) +! +! integer :: olog = 10 +! +! real(kind=DP) :: ring_a = 0.5_DP + + real(kind=DP) :: ring_a + + public ring_a + + +END MODULE ring_header + + +MODULE ring_func + + use ring_header + use GKV_math, only: math_eli1, math_eli2 + + implicit none + + private + + public func_k, func_g, func_psi, func_x, func_eli1, func_eli2 + + +CONTAINS + + + FUNCTION func_k( r, z ) + + real(kind=DP) :: func_k + real(kind=DP), intent(in) :: r, z + + func_k = sqrt( ( 4._DP * ring_a * abs(r) ) / ( ( ring_a + r )**2 + z**2 ) ) + + END FUNCTION func_k + + + FUNCTION func_g( r, z ) + + real(kind=DP) :: func_g + real(kind=DP), intent(in) :: r, z + + real(kind=DP) :: eli1, eli2, kk + + kk = func_k( r, z )**2 + + call math_eli1( kk, eli1 ) +!!! if( icon /= 0 ) print *, "# icon in celi1 = ", icon + call math_eli2( kk, eli2 ) +!!! if( icon /= 0 ) print *, "# icon in celi2 = ", icon + + func_g = ( ( 1._DP - 0.5_DP * kk ) * eli1 - eli2 ) & + * sqrt( ( ring_a + r )**2 + z**2 ) * 0.5_DP + + END FUNCTION func_g + + + FUNCTION func_psi( r, z ) + + real(kind=DP) :: func_psi + real(kind=DP), intent(in) :: r, z + + func_psi = func_g( r, z ) / func_g( 1._DP, 0._DP ) + + END FUNCTION func_psi + + + FUNCTION func_x( r, z, psi0 ) + + real(kind=DP) :: func_x + real(kind=DP), intent(in) :: r, z, psi0 + + func_x = (func_psi( 1._DP, 0._DP ) - func_psi( r, z ) ) / psi0 + + END FUNCTION func_x + + + FUNCTION func_eli1( r, z ) + + real(kind=DP) :: func_eli1 + real(kind=DP), intent(in) :: r, z + + real(kind=DP) :: eli1, kk + + kk = func_k( r, z )**2 + + call math_eli1( kk, eli1 ) + + func_eli1 = eli1 + + END FUNCTION func_eli1 + + + FUNCTION func_eli2( r, z ) + + real(kind=DP) :: func_eli2 + real(kind=DP), intent(in) :: r, z + + real(kind=DP) :: eli2, kk + + kk = func_k( r, z )**2 + + call math_eli2( kk, eli2 ) + + func_eli2 = eli2 + + END FUNCTION func_eli2 + + +END MODULE ring_func + + +MODULE ring_diff + + use ring_header + + implicit none + + private + + public diff_r, diff_z, diff_rho + + +CONTAINS + + + FUNCTION diff_r( fun, rin, zin ) + + real(kind=DP), external :: fun + real(kind=DP), intent(in) :: rin, zin + + real(kind=DP) :: diff_r + + real(kind=DP) :: dr1, dr2 + + + dr1 = abs( rin ) * 1.d-4 + dr2 = dr1 * 2._DP + + if( rin == 0._DP ) then + diff_r = 0._DP + + else + diff_r = ( - fun(rin+dr2,zin) + 8._DP*fun(rin+dr1,zin) & + - 8._DP*fun(rin-dr1,zin) + fun(rin-dr2,zin) ) & + / ( 12._DP * dr1 ) + + end if + + + END FUNCTION diff_r + + + FUNCTION diff_z( fun, rin, zin ) + + real(kind=DP), external :: fun + real(kind=DP), intent(in) :: rin, zin + + real(kind=DP) :: diff_z + + real(kind=DP) :: dz1, dz2 + + + dz1 = abs( zin ) * 1.d-4 + dz2 = dz1 * 2._DP + +! if( zin == 0._DP ) then + if( abs(zin) < 1.d-12) then + diff_z = 0._DP + + else + diff_z = ( - fun(rin,zin+dz2) + 8._DP*fun(rin,zin+dz1) & + - 8._DP*fun(rin,zin-dz1) + fun(rin,zin-dz2) ) & + / ( 12._DP * dz1 ) + + end if + + + END FUNCTION diff_z + + + FUNCTION diff_rho( fun, rin, zin ) + + real(kind=DP), external :: fun + real(kind=DP), intent(in) :: rin, zin + + real(kind=DP) :: diff_rho + +! real(kind=DP) :: rho, tht + real(kind=DP) :: tht + + +! rho = sqrt( ( rin - ring_a )**2 + zin**2 ) + tht = atan2( zin, rin - ring_a ) + + diff_rho = diff_r( fun, rin, zin ) * cos( tht ) & + + diff_z( fun, rin, zin ) * sin( tht ) + + + END FUNCTION diff_rho + + +END MODULE ring_diff + + +MODULE ring_bfld + + use ring_header + use ring_func + use ring_diff + + implicit none + + private + + public bfld_br, bfld_bz, bfld_magb, bfld_gradbr, bfld_gradbz + + +CONTAINS + + + FUNCTION bfld_br( r, z ) + + real(kind=DP) :: bfld_br + real(kind=DP), intent(in) :: r, z + +! bfld_br = diff_z( func_psi, r, z) / r + + bfld_br = diff_z( func_psi, r, z ) / r + + END FUNCTION bfld_br + + + FUNCTION bfld_bz( r, z ) + + real(kind=DP) :: bfld_bz + real(kind=DP), intent(in) :: r, z + +! bfld_bz = diff_r( func_psi, r, z ) / r + + bfld_bz = - diff_r( func_psi, r, z ) / r + + END FUNCTION bfld_bz + + + FUNCTION bfld_magb( r, z ) + + real(kind=DP) :: bfld_magb + real(kind=DP), intent(in) :: r, z + + real(kind=DP) :: br, bz + + br = bfld_br( r, z ) + bz = bfld_bz( r, z ) + + bfld_magb = sqrt( br**2 + bz**2 ) + +!!! print *, br, bz + + END FUNCTION bfld_magb + + + FUNCTION bfld_gradbr( r, z ) + + real(kind=DP) :: bfld_gradbr + real(kind=DP), intent(in) :: r, z + + bfld_gradbr = diff_r( bfld_magb, r, z ) + + END FUNCTION bfld_gradbr + + + FUNCTION bfld_gradbz( r, z ) + + real(kind=DP) :: bfld_gradbz + real(kind=DP), intent(in) :: r, z + + bfld_gradbz = diff_z( bfld_magb, r, z ) + + END FUNCTION bfld_gradbz + + +END MODULE ring_bfld + + +MODULE GKV_ring +!----------------------------------------------------------------- +! +! Flux tube coordinates in the ring dipole geometry +! +! Definition of the flux tube coordinates +! in the ring dipole geometry +! +! x = (Psi_0 - Psi) / (R_0*B'_0) +! y = R_0 * phi +! z = Theta (= arctan(Z/(R-a)) +! where (R, phi, Z ) are the cylindorical coordinates +! +! B'_0 = 1/R * (dPsi/dr - dPsi/dz) | r=1,z=0 +! B_0 = B / B'_0 | r=1,z=0 +! We also use the normalization of Psi_0*(gradx crs grady) = B_0*R_0^2 +! with R_0 = 1 and B_0 = 1 as the units +! +!----------------------------------------------------------------- + + use ring_header +! use gkv_header + use ring_func + use ring_diff + use ring_bfld + + implicit none + + private + + public ring_coordinates + + +CONTAINS + + SUBROUTINE ring_coordinates( a, tht, bb, ub_dot_grdb, ub_crs_grdb, & + gxx, gxy, gxz, gyy, gyz, gzz, rootg, dbdx, dbdz) + + real(kind=DP), intent(in) :: a, tht + real(kind=DP), intent(out) :: bb, ub_dot_grdb, ub_crs_grdb + real(kind=DP), intent(out) :: gxx, gxy, gxz, gyy, gyz, gzz, rootg +!>> + real(kind=DP), intent(out) :: dbdx, dbdz + real(kind=DP) :: R0, psi0 + real(kind=DP) :: eps_x, rho_p, rho_m, r_p, r_m, z_p, z_m +!<< + real(kind=DP) :: r, z, rho, rho1, hh, dh + + real(kind=DP) :: b0 + real(kind=DP) :: gbr, gbz, ubr, ubz!, psi_n, b_rootg_i, ub_dot_grdh + +! real(kind=DP) :: rootg2 + + real(kind=DP) :: eps, eps0 = 0.00000001_DP + integer :: ic, nc = 20 + + + ring_a = a + + rho = 1._DP - ring_a + + hh = 0._DP + dh = acos(-1._DP) * 0.01_DP + + + +! compute for hh +!!! do while ( hh < abs(tht)-dh*0.5 ) + do while ( hh < abs(tht)-dh ) + hh = hh + dh + + eps = 1._DP + ic = 0 + + r = rho * cos(hh) + ring_a + z = rho * sin(hh) + + do while ( eps > eps0 .AND. ic < nc ) + + ic = ic + 1 + + rho1 = rho - ( func_psi( r, z ) - 1._DP ) / diff_rho( func_psi, r, z ) + eps = abs( rho1 - rho ) / abs( rho ) + rho = rho1 + + r = rho * cos(hh) + ring_a + z = rho * sin(hh) + if( abs(z) < 1.d-14 ) z = 0._DP + + end do +!!! print *, "# ic, tht, rho1, psi = ", ic, tht, rho, func_psi( r, z ) + + end do + +! compute for tht + eps = 1._DP + ic = 0 + + r = rho * cos(tht) + ring_a + z = rho * sin(tht) + + do while ( eps > eps0 .AND. ic < nc ) + + ic = ic + 1 + + rho1 = rho - ( func_psi( r, z ) - 1._DP ) / diff_rho( func_psi, r, z ) + eps = abs( rho1 - rho ) / abs( rho ) + rho = rho1 + + r = rho * cos(tht) + ring_a + z = rho * sin(tht) + if( abs(z) < 1.d-14 ) z = 0._DP + + end do + + if( ic == nc ) then + print *, "# ic, tht, rho1, psi = ", ic, tht, rho, func_psi( r, z ) + end if + +!>> + R0 = 1._DP + eps_x = 0.0000001_DP + rho_p = rho + eps_x + r_p = rho_p * cos(tht) + ring_a + z_p = rho_p * sin(tht) + rho_m = rho - eps_x + r_m = rho_m * cos(tht) + ring_a + z_m = rho_m * sin(tht) +!<< + + b0 = bfld_magb( 1._DP, 0._DP ) + + bb = bfld_magb ( r, z )/b0 + +!>> + psi0 = b0*R0**2 + dbdx = ( bfld_magb( r_p, z_p ) - bfld_magb( r_m, z_m ) )/b0/( func_x(r_p, z_p, psi0) - func_x(r_m, z_m, psi0) ) +!<< + + ubr = bfld_br( r, z )/b0 / bb + ubz = bfld_bz( r, z )/b0 / bb + + gbr = bfld_gradbr( r, z )/b0 + gbz = bfld_gradbz( r, z )/b0 + + ub_dot_grdb = ubr * gbr + ubz * gbz + ub_crs_grdb = ubz * gbr - ubr * gbz + +! ub_dot_grdh = - ubr * sin( tht ) + ubz * cos( tht ) + +! psi_n = func_psi( r, z ) / func_psi( 1._DP, 0._DP ) + + gxx = ( r * bb )**2 + gxy = 0._DP + gxz = - r * bb * ( ubz*sin(tht) + ubr*cos(tht) ) / rho + gyy = 1._DP / r**2 + gyz = 0._DP + gzz = 1._DP / rho**2 + rootg = 1._DP / sqrt( gyy*( gxx*gzz - gxz**2 ) ) + +!>> + dbdz = ub_dot_grdb * bb * rootg +!<< + +! b_rootg_i = 1._DP / ( bb * rootg ) + +! rootg2= r * rho * b0 / ( diff_r( func_psi, r, z )*cos(tht) & +! + diff_z( func_psi, r, z )*sin(tht) ) +! rootg2= 1._DP / ( ubz * bb * cos(tht) & +! - ubr * bb * sin(tht) ) + +!! debug +! write(unit=6, fmt="(1p, 32e15.7)" ) tht, rho, r, z, func_psi(r,z), & +! bb, ubr, ubz, gbr, gbz, ub_dot_grdb, ub_crs_grdb, ub_dot_grdh, & +! psi_n, gxx, gxy, gxz, gyy, gyz, gzz, rootg, b_rootg_i, rootg2 +!! debug + + + + END SUBROUTINE ring_coordinates + + +END MODULE GKV_ring diff --git a/src/gkvp_set.f90 b/src/gkvp_set.f90 index 37259b9..c541732 100644 --- a/src/gkvp_set.f90 +++ b/src/gkvp_set.f90 @@ -5,6 +5,9 @@ MODULE GKV_set ! ! Update history of gkvp_set.f90 ! -------------- +! gkvp_f0.62 (S. Maeyama, Mar 2023) +! - Contents of subroutine set_cnfig are moved to GKV_geom, to implement +! time-dependent metrics and operators in rotating flux-tube model. ! gkvp_f0.61 (S. Maeyama, Mar 2021) ! - equib_type = "s-alpha-shift" is added. ! - Initial random phase rr is set by global (mx,gmy) indices. @@ -21,24 +24,20 @@ MODULE GKV_set use GKV_header use GKV_mpienv - use GKV_math, only: math_j0, math_j1, math_j2, math_g0, math_random - use GKV_intgrl, only: intgrl_fsrf, intgrl_v0_moment_ms + use GKV_math, only: math_random use GKV_fld, only: fld_esfield, fld_emfield_ff, fld_ff2hh use GKV_bndry, only: bndry_zvm_bound_f use GKV_advnc, only: caldlt_rev use GKV_dtc, only: dtc_init -! for vmec equilibrium -! use GKV_vmecin, only: vmecin_fileopen, vmecin_coeff, vmecin_read -! for vmec equilibrium w/ Booz_xform by M. Nakata & M. Nunami (Aug. 2016) - use GKV_vmecbzx, only: vmecbzx_boozx_read, vmecbzx_boozx_coeff -! for tokamak(eqdsk) equilibrium - use GKV_igs, only: igs_read, igs_coeff use GKV_colli, only: colli_set_param use GKV_colliimp, only: colliimp_set_param use GKV_tips, only: tips_reality !fj start 202010 use GKV_fileio !fj end 202010 + use GKV_geom, only : geom_read_nml, geom_init_kxkyzvm, & + geom_init_metric, geom_set_operators, & + geom_reset_time implicit none @@ -46,7 +45,6 @@ MODULE GKV_set public set_init, set_close - CONTAINS !-------------------------------------- @@ -72,7 +70,8 @@ SUBROUTINE set_init( ff, phi, Al, hh, time ) trim(equib_type) == "s-alpha-shift" .OR. & trim(equib_type) == "circ-MHD" .OR. & trim(equib_type) == "vmec" .OR. & - trim(equib_type) == "eqdsk" ) then + trim(equib_type) == "eqdsk" .OR. & + trim(equib_type) == "ring" ) then call set_cnfig @@ -185,6 +184,7 @@ SUBROUTINE set_start if( rankg == 0 ) then open( omtr, file=trim(f_hst)//"mtr."//cnew ) + open( omtf, file=trim(f_hst)//"mtf."//cnew ) open( odtc, file=trim(f_hst)//"dtc."//cnew ) open( oeng, file=trim(f_hst)//"eng."//cnew ) open( omen, file=trim(f_hst)//"men."//cnew ) @@ -271,6 +271,7 @@ SUBROUTINE set_close if( rankg == 0 ) then close( omtr ) + close( omtf ) close( odtc ) close( oeng ) close( omen ) @@ -394,99 +395,10 @@ END SUBROUTINE set_param SUBROUTINE set_cnfig !-------------------------------------- - real(kind=DP) :: s_hat - - real(kind=DP) :: eps_r - - real(kind=DP) :: rdeps00, eps_hor, lprd, mprd, lmmq, malpha - real(kind=DP) :: eps_mor, eps_por, lprdm1, lprdp1, lmmqm1, lmmqp1 - real(kind=DP) :: eps_rnew, rdeps1_0, rdeps1_10, rdeps2_10, rdeps3_10 - -! for s-alpha model with Shafranov shift - real(kind=DP) :: p_total, dp_totaldx, beta_total, alpha_MHD - -! for circular MHD - real(kind=DP), dimension(1:3,1:3) :: gg - real(kind=DP) :: kkx, kky, domgdz, domgdx, domgdy - - -!! for vmec equilibrium -! real(kind=DP) :: rho2R_0, q_input, theta -! real(kind=DP) :: r_0 -! real(kind=DP) :: gdwss, gdwtt, gdwzz, gdwst, gdwsz, gdwtz, & -! gupss, guptt, gupzz, gupst, gupsz, guptz, & -! babs, Bs , Bth , Bzt , dBds, dBdt, dBdz, & -! dBdt_mir, vmec_rootg, rootgft, rootgbz - real(kind=DP) :: theta - - - real(kind=DP) :: lx, ly, lz, kxmin, kymin, dz, mmax, dm, del_c - real(kind=DP) :: lz_l, z0, z0_l - integer :: n_tht, m_j - - real(kind=DP) :: gg0 - - real(kind=DP) :: bb, kmo - real(kind=DP) :: cfsrf_l - - integer :: global_iv, global_im - integer :: mx, my, iz, iv, im, is, is1, is2, ierr_mpi - - - complex(kind=DP), dimension(:,:,:,:,:), allocatable :: wf - complex(kind=DP), dimension(:,:,:), allocatable :: nw - real(kind=DP), dimension(:,:,:), allocatable :: ww - -! real(kind=DP) :: rad_a, r_minor, eps_b, rho_unit, r_a -! real(kind=DP) :: R0_unit, r_edge, b0b00, alpha_fix - - real(kind=DP), dimension(0:ns-1) :: eta real(kind=DP), dimension(0:ns-1,0:ns-1) :: nust - real(kind=DP) :: r_major + real(kind=DP) :: lx, ly, eps_r + integer :: is1, is2 - real(kind=DP) :: s_input, s_0 ! radial label of fluxtube center - integer :: mc_type ! 0:Axisym., 1:Boozer, 2:Hamada - integer :: q_type ! 0:use q and s_hat value in confp, 1:calclated by IGS - integer :: isw, nss, ntheta, nzeta - real(kind=DP) :: phi_ax ! axisymetric toroidal angle - - integer, parameter :: num_omtr = 13 - real(kind=DP) :: metric_l(1:num_omtr,-nz:nz-1), metric_g(1:num_omtr,-global_nz:global_nz-1) - - - - - namelist /physp/ R0_Ln, & ! R0/Lns - R0_Lt, & ! R0/Lts - nu, & ! factor for collision freq. in LB model - Anum, & ! mass number - Znum, & ! charge number - fcs, & ! charge-density fraction - sgn, & ! signs of charge - tau, & ! T-ratio Ts/T0, T0=reference ion temp. of ranks=1 - dns1, & ! initial perturbation amplitude - tau_ad, & ! Ti/Te for ITG-ae, Te/Ti for ETG-ai - lambda_i, & ! (Debye/rho_tp)^2 - beta, & ! mu0*ni*Ti/B^2 - ibprime, & ! flag for finite beta-prime effect on kvd - vmax, & ! maximum v_para in unit of v_ts - nx0 ! mode number for the initial perturbation - - namelist /rotat/ mach, uprime, gamma_e - - namelist /nperi/ n_tht, kymin, m_j, del_c - namelist /confp/ eps_r, eps_rnew, & - q_0, s_hat, & - lprd, mprd, eps_hor, eps_mor, eps_por, & - rdeps00, rdeps1_0, rdeps1_10, & - rdeps2_10, rdeps3_10, malpha -! namelist /vmecp/ q_0, rad_a, & -! R0_unit, r_edge, & -! b0b00, alpha_fix - namelist /vmecp/ s_input, nss, ntheta, nzeta - - namelist /igsp/ s_input, mc_type, q_type, nss, ntheta - namelist /nu_ref/ Nref, & ! reference (electron) density in m^(-3) Lref, & ! reference length (=R_axis) in m Tref, & ! reference main-ion (ranks=1) temperature in keV @@ -494,1309 +406,15 @@ SUBROUTINE set_cnfig iFLR, & ! flag for GK- or DK-limit in collision icheck ! flag for Maxwellain anihilation test (w/ iFLR=0) - tau(:) = 1.0_DP - nu(:) = 0.002_DP - R0_Ln(:) = 2.5_DP - R0_Lt(:) = 7.5_DP - - - read(inml,nml=physp) - - - do is = 0, ns-1 - if( R0_Ln(is) /= 0._DP ) then - eta(is) = R0_Lt(is) / R0_Ln(is) - else - eta(is) = 1.d+20 - end if - end do - - - write( olog, * ) " # Physical parameters" - write( olog, * ) "" - write( olog, * ) " # r_major/L_ns = ", R0_Ln(:) - write( olog, * ) " # r_major/L_ts = ", R0_Lt(:) - write( olog, * ) " # eta = ", eta(:) - write( olog, * ) " # nu = ", nu(:) - write( olog, * ) " # A-number = ", Anum(:) - write( olog, * ) " # Z-number = ", Znum(:) - write( olog, * ) " # fcs = ", fcs(:) - write( olog, * ) " # sgn = ", sgn(:) - write( olog, * ) " # tau = ", tau(:) - write( olog, * ) " # dns1 = ", dns1(:) - write( olog, * ) " # tau_ad = ", tau_ad - write( olog, * ) " # lambda_i^2 = ", lambda_i - write( olog, * ) " # beta_i = ", beta - write( olog, * ) " # ibprime = ", ibprime - write( olog, * ) " # nx0 = ", nx0 - write( olog, * ) "" - - - mach = 0._DP - uprime = 0._DP - gamma_e = 0._DP - - read(inml,nml=rotat) - - write( olog, * ) " # Mean rotation parameters" - write( olog, * ) "" - write( olog, * ) " # Mach number = ", mach - write( olog, * ) " # uptime = ", uprime - write( olog, * ) " # gamma_ExB = ", gamma_e - write( olog, * ) "" - - - n_tht = 1 - - read(inml,nml=nperi) - - - if( trim(equib_type) == "slab") then - - read(inml,nml=confp) - - lprdm1 = 0._DP - lprdp1 = 0._DP - - lmmq = 0._DP - lmmqm1 = 0._DP - lmmqp1 = 0._DP - - q_0 = 1._DP ! For now, fixed q_0=1. Changing q_0 can extend parallel z-box size. - s_hat = 0._DP ! only shear less slab - eps_r = 1._DP - - eps_hor = 0._DP - lprd = 0._DP - mprd = 0._DP - malpha = 0._DP - - rdeps00 = 0._DP - eps_mor = 0._DP - eps_por = 0._DP - - write( olog, * ) " # Configuration parameters" - write( olog, * ) "" - write( olog, * ) " # q_0 = ", q_0 - write( olog, * ) " # s_hat = ", s_hat - write( olog, * ) " # eps_r = ", eps_r - write( olog, * ) "" - - write( olog, * ) " # eps_hor = ", eps_hor - write( olog, * ) " # lprd = ", lprd - write( olog, * ) " # mprd = ", mprd - write( olog, * ) " # malpha = ", malpha - write( olog, * ) " # rdeps00 = ", rdeps00 - - write( olog, * ) " # eps_mor = ", eps_mor - write( olog, * ) " # lprdm1 = ", lprdm1 - write( olog, * ) " # eps_por = ", eps_por - write( olog, * ) " # lprdp1 = ", lprdp1 - write( olog, * ) "" - - else if( trim(equib_type) == "analytic" .OR. & - trim(equib_type) == "s-alpha" .OR. & - trim(equib_type) == "s-alpha-shift" .OR. & - trim(equib_type) == "circ-MHD" ) then - - - read(inml,nml=confp) - - - lprdm1 = lprd - 1.0_DP - lprdp1 = lprd + 1.0_DP - - lmmq = lprd - mprd * q_0 - lmmqm1 = lprdm1 - mprd * q_0 - lmmqp1 = lprdp1 - mprd * q_0 - - - write( olog, * ) " # Configuration parameters" - write( olog, * ) "" - write( olog, * ) " # q_0 = ", q_0 - write( olog, * ) " # s_hat = ", s_hat - write( olog, * ) " # eps_r = ", eps_r - write( olog, * ) "" - - write( olog, * ) " # eps_hor = ", eps_hor - write( olog, * ) " # lprd = ", lprd - write( olog, * ) " # mprd = ", mprd - write( olog, * ) " # malpha = ", malpha - write( olog, * ) " # rdeps00 = ", rdeps00 - - write( olog, * ) " # eps_mor = ", eps_mor - write( olog, * ) " # lprdm1 = ", lprdm1 - write( olog, * ) " # eps_por = ", eps_por - write( olog, * ) " # lprdp1 = ", lprdp1 - write( olog, * ) "" - - - -! else if( trim(equib_type) == "vmec" ) then -! -! -!! --- Paramters at rho=0.65 (shot#088343 at t = 1.833 [s]) -!! -!! Ln_unit =-4.230701_DP ! Ln [m] -!! Lt_unit = 0.3135611_DP ! Lt [m] -!! R0_unit = 3.599858_DP ! R0 [m] -!! r_edge = 0.6362872D0 ! r_edge [m] -!! b0b00 = 2.88846853946973647d0/2.940307D0 ! b00mode/B0 -!! alpha_fix = 0.314159253589793d0 ! pi/10 -!! -! -! read(inml,nml=vmecp) -! -! eps_b = r_edge / R0_unit ! --- a / R ! by nunami (2010.04.21) -! -! rho2R_0 = eps_b / rad_a ! --- rho / R_0 -! rho_unit = rho2R_0 * R0_unit ! --- rho -! r_a = rad_a * rho2R_0 ! --- rad_a / R_0 -! -! eps_r = 0.1115200537d0 -! -! call vmecin_fileopen -! -! -! call vmecin_read -! -! -! q_input = q_0 -! theta = 0._DP -! -! call vmecin_coeff( rad_a, R0_unit, rho2R_0, q_input, theta, & -! alpha_fix, r_0, r_minor, s_hat, & -! gdwss, gdwtt, gdwzz, gdwst, gdwsz, gdwtz, & -! gupss, guptt, gupzz, gupst, gupsz, guptz, & -! babs, Bs, Bth, Bzt, dBds, dBdt, dBdz, & -! dBdt_mir, vmec_rootg, rootgft, rootgbz ) -! -! -! write( olog, * ) " # Configuration parameters" -! write( olog, * ) "" -! write( olog, * ) " # q_0 = ", q_0 -! write( olog, * ) " # s_hat = ", s_hat -! write( olog, * ) "" -! write( olog, * ) " # eps_r = ", eps_r -! write( olog, * ) "" - - - else if( trim(equib_type) == "vmec" ) then - - - read(inml,nml=confp) - - read(inml,nml=vmecp) - - call vmecbzx_boozx_read( nss, ntheta, nzeta ) - - isw = 0 - iz = 0 - call vmecbzx_boozx_coeff( isw, nss, ntheta, nzeta, s_input, iz, 0._DP, lz_l, & ! input - s_0, q_0, s_hat, eps_r, phi_ax, & ! output - omg(iz), rootg(iz), domgdx, domgdz, domgdy, & - gg(1,1), gg(1,2), gg(1,3), gg(2,2), & - gg(2,3), gg(3,3) ) - - - - write( olog, * ) " # Configuration parameters" - write( olog, * ) "" - write( olog, * ) " # r_major/L_ns = ", R0_Ln(:) - write( olog, * ) " # r_major/L_ts = ", R0_Lt(:) - write( olog, * ) " # eta = ", eta(:) - write( olog, * ) " # q_0 = ", q_0 - write( olog, * ) " # s_hat = ", s_hat - write( olog, * ) " # eps_r = ", eps_r - write( olog, * ) " # s_input, s_0 = ", s_input, s_0 - write( olog, * ) " # nss, ntheta, nzeta = ", nss, ntheta, nzeta - - - else if( trim(equib_type) == "eqdsk" ) then - - - read(inml,nml=confp) - - read(inml,nml=igsp) - - call igs_read( mc_type, nss, ntheta ) - - if ( q_type == 1 ) then - isw = 0 - iz = 0 - call igs_coeff( isw, mc_type, nss, ntheta, s_input, 0._DP, lz_l, & ! input - s_0, q_0, s_hat, eps_r, theta, & ! output - omg(iz), rootg(iz), domgdx, domgdz, domgdy, & - gg(1,1), gg(1,2), gg(1,3), gg(2,2), & - gg(2,3), gg(3,3) ) - end if - - - - write( olog, * ) " # Configuration parameters" - write( olog, * ) "" - write( olog, * ) " # r_major/L_ns = ", R0_Ln(:) - write( olog, * ) " # r_major/L_ts = ", R0_Lt(:) - write( olog, * ) " # eta = ", eta(:) - write( olog, * ) " # q_0 = ", q_0 - write( olog, * ) " # s_hat = ", s_hat - write( olog, * ) " # eps_r = ", eps_r - write( olog, * ) " # s_input, s_0 = ", s_input, s_0 - write( olog, * ) " # nss, ntheta = ", nss, ntheta - - else - - write( olog, * ) " # wrong choice of the equilibrium " - call flush(olog) - call MPI_Finalize(ierr_mpi) - stop - - end if - - -! --- coordinate settings --- - - if (abs(s_hat) < 1.d-10) then ! When s_hat == ZERO - m_j = 0 - kxmin = kymin - else if (m_j == 0) then - kxmin = kymin - else - kxmin = abs(2._DP * pi * s_hat * kymin / real(m_j, kind=DP)) - end if - lx = pi / kxmin - ly = pi / kymin - ! kymin=pi/ly=pi/[r_minor*pi/(q0*n_alp)]=q0*n_alp/r_minor - - lz = real( n_tht, kind=DP ) * pi ! total z-length - lz_l = lz / real( nprocz, kind=DP ) ! local z-length - - do mx = -nx, nx - kx(mx) = kxmin * real( mx, kind=DP ) - end do - - ky(:) = 0._DP - do my = ist_y_g, iend_y_g - ky(my-ist_y_g) = kymin * real( my, kind=DP ) - end do - - kxmin_g = kxmin - kymin_g = kymin - - z0 = - lz ! global lower boundary - z0_l = 2._DP * lz_l * real( rankz, kind=DP ) + z0 - ! local lower boundary - - dz = lz_l / real( nz, kind=DP ) - - do iz = -nz, nz-1 - zz(iz) = dz * real( iz + nz, kind=DP ) + z0_l - end do - - - dv = 2._DP * vmax / real( 2 * nv * nprocv -1, kind=DP ) - - do iv = 1, 2*nv - global_iv = 2 * nv * rankv + iv - vl(iv) = dv * ( real( global_iv - nv * nprocv - 1, kind=DP ) + 0.5_DP ) - end do - ! --- debug - ! write( olog, * ) " *** iv, vl " - ! do iv = 1, 2*nv - ! global_iv = 2 * nv * rankv + iv - ! write( olog, * ) iv, global_iv, vl(iv) - ! end do - ! write( olog, * ) "" - - mmax = vmax - 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 - - - write( olog, * ) " # Numerical parameters" - write( olog, * ) "" - write( olog, * ) " # n_tht = ", n_tht - write( olog, * ) " # lx, ly, lz = ", lx, ly, lz - write( olog, * ) " # lz, z0 = ", lz, z0 - write( olog, * ) " # lz_l, z0_l = ", lz_l, z0_l - write( olog, * ) " # kxmin, kymin = ", kxmin, kymin - write( olog, * ) " # kxmax, kymax = ", kxmin*nx, kymin*global_ny - write( olog, * ) " # kperp_max = ", sqrt((kxmin*nx)**2+(kymin*global_ny)**2) - write( olog, * ) " # m_j, del_c = ", m_j, del_c - write( olog, * ) " # dz = ", dz - write( olog, * ) " # dv, vmax = ", dv, vmax - write( olog, * ) " # dm, mmax = ", dm, mmax - write( olog, * ) "" - - if (gamma_e == 0._DP) then - tlim_exb = 999999.d0 - else - tlim_exb = (kxmin*(nx-nx0))/(kymin*global_ny*abs(gamma_e)) - end if - write( olog, * ) " # ExB limit time tlim_exb = ", tlim_exb - write( olog, * ) " # for (mx=nx0,my=global_ny) initial perturbation: " - write( olog, * ) " # tlim_exb = kxmin*(nx-nx0)/(kymax*|gamma_e|)" - write( olog, * ) "" - - -! --- coordinate settings --- - - -! --- operator settings --- - - - do iz = -nz, nz-1 - -!!! for slab model - if ( trim(equib_type) == "slab") then - - q_bar = q_0 - r_major = 1._DP ! in the R0 unit - theta = zz(iz) - - omg(iz) = 1._DP - rootg(iz) = q_0*r_major - dpara(iz) = dz * q_0 * r_major - - do im = 0, nm - vp(iz,im) = sqrt( 2._DP * mu(im) )!* omg(iz) ) - mir(iz,im) = 0._DP - do iv = 1, 2*nv - vdx(iz,iv,im) = 0._DP - vdy(iz,iv,im) = 0._DP - 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 ) ) - end do - end do ! im loop ends - - ksq(:,:,iz) = 0._DP - do my = ist_y, iend_y - do mx = -nx, nx - ksq(mx,my,iz) = kx(mx)**2 + ky(my)**2 - end do - end do - - baxfactor = 1._DP - - !- for OUTPUT hst/*.mtr.* - - domgdz = 0._DP - domgdy = 0._DP - domgdx = 0._DP - gg(1,1) = 1._DP - gg(1,2) = 0._DP - gg(1,3) = 0._DP - gg(2,1) = gg(1,2) - gg(2,2) = 1._DP - gg(2,3) = 0._DP - gg(3,1) = gg(1,3) - gg(3,2) = gg(2,3) - gg(3,3) = 1._DP - metric_l( 1,iz) = zz(iz) ! [ 1] - metric_l( 2,iz) = theta ! [ 2] - metric_l( 3,iz) = omg(iz) ! [ 3] - metric_l( 4,iz) = domgdx ! [ 4] - metric_l( 5,iz) = domgdy ! [ 5] - metric_l( 6,iz) = domgdz ! [ 6] - metric_l( 7,iz) = gg(1,1) ! [ 7] - metric_l( 8,iz) = gg(1,2) ! [ 8] - metric_l( 9,iz) = gg(1,3) ! [ 9] - metric_l(10,iz) = gg(2,2) ! [10] - metric_l(11,iz) = gg(2,3) ! [11] - metric_l(12,iz) = gg(3,3) ! [12] - metric_l(13,iz) = rootg(iz)! [13] - !------------------------- - - -!!! for the concentric and large-aspect-ratio model !!! - else if( trim(equib_type) == "analytic" ) then - - q_bar = q_0 - r_major = 1._DP ! in the R0 unit - - theta = zz(iz) - - omg(iz) = 1._DP & - - eps_r * ( cos( zz(iz) ) & - + eps_hor * cos( lmmq * zz(iz) - malpha ) & - + eps_mor * cos( lmmqm1 * zz(iz) - malpha ) & - + eps_por * cos( lmmqp1 * zz(iz) - malpha ) ) - - rootg(iz) = q_0*r_major/omg(iz) - dpara(iz) = dz * q_0 * r_major - - ! --- debug - ! write( olog, * ) " *** z, omg " - ! do iz = -nz, nz-1 - ! write( olog, * ) zz(iz), omg(iz) - ! end do - ! write( olog, * ) - - 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 ) & - * ( sin(zz(iz)) & - + eps_hor * lmmq * sin( lmmq * zz(iz) - malpha ) & - + eps_mor * lmmqm1 * sin( lmmqm1 * zz(iz) - malpha ) & - + eps_por * lmmqp1 * sin( lmmqp1 * zz(iz) - malpha ) ) - - do iv = 1, 2*nv - !do my = ist_y, iend_y - ! do mx = -nx, nx - ! kvd and kvs are revised November 2011 - ! into general species forms. - !!!!kvd(mx,my,iz,iv,im)= & - !!!! - ( vl(iv)**2 + omg(iz)*mu(im) ) * eps_rnew / r_major & - !!!! * ( ky(my) * ( rdeps00 + rdeps1_0 * cos( zz(iz) ) & - !!!! + rdeps2_10 * cos( lmmq * zz(iz) - malpha ) & - !!!! + rdeps1_10 * cos( lmmqm1 * zz(iz) - malpha ) & - !!!! + rdeps3_10 * cos( lmmqp1 * zz(iz) - malpha ) ) & - !!!! + ( kx(mx) + s_hat * zz(iz) * ky(my) ) & - !!!! * ( sin( zz(iz) ) & - !!!! + eps_hor * lprd * sin( lmmq * zz(iz) - malpha ) & - !!!! + eps_mor * lprdm1 * sin( lmmqm1 * zz(iz) - malpha ) & - !!!! + eps_por * lprdp1 * sin( lmmqp1 * zz(iz) - malpha ) ) & - !!!! ) * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) - - !!!!kvs(mx,my,iz,iv,im) = & - !!!! - sgn(ranks) * tau(ranks) / Znum(ranks) * ky(my) & - !!!! * ( R0_Ln(ranks) + R0_Lt(ranks) * ( 0.5_DP*vl(iv)**2 & - !!!! + omg(iz)*mu(im) - 1.5_DP ) ) - - !%%% kvd = kx(mx) * vdx(iz,iv,im) + ky(my) * vdy(iz,iv,im) %%% - !%%% kvs = ky(my) * vsy(iz,iv,im) %%% - vdx(iz,iv,im)= & - - ( vl(iv)**2 + omg(iz)*mu(im) ) * eps_rnew / r_major & - * ( 0._DP * ( rdeps00 + rdeps1_0 * cos( zz(iz) ) & - + rdeps2_10 * cos( lmmq * zz(iz) - malpha ) & - + rdeps1_10 * cos( lmmqm1 * zz(iz) - malpha ) & - + rdeps3_10 * cos( lmmqp1 * zz(iz) - malpha ) ) & - + ( 1._DP + s_hat * zz(iz) * 0._DP ) & - * ( sin( zz(iz) ) & - + eps_hor * lprd * sin( lmmq * zz(iz) - malpha ) & - + eps_mor * lprdm1 * sin( lmmqm1 * zz(iz) - malpha ) & - + eps_por * lprdp1 * sin( lmmqp1 * zz(iz) - malpha ) ) & - ) * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) - - vdy(iz,iv,im)= & - - ( vl(iv)**2 + omg(iz)*mu(im) ) * eps_rnew / r_major & - * ( 1._DP * ( rdeps00 + rdeps1_0 * cos( zz(iz) ) & - + rdeps2_10 * cos( lmmq * zz(iz) - malpha ) & - + rdeps1_10 * cos( lmmqm1 * zz(iz) - malpha ) & - + rdeps3_10 * cos( lmmqp1 * zz(iz) - malpha ) ) & - + ( 0._DP + s_hat * zz(iz) * 1._DP ) & - * ( sin( zz(iz) ) & - + eps_hor * lprd * sin( lmmq * zz(iz) - malpha ) & - + eps_mor * lprdm1 * sin( lmmqm1 * zz(iz) - malpha ) & - + eps_por * lprdp1 * sin( lmmqp1 * zz(iz) - malpha ) ) & - ) * ( 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 ) ) - - ! end do - !end do - end do - - end do ! im loop ends - - - ksq(:,:,iz) = 0._DP - do my = ist_y, iend_y - do mx = -nx, nx - ksq(mx,my,iz) = ( kx(mx) + s_hat * zz(iz) * ky(my) )**2 + ky(my)**2 - end do - end do - - baxfactor = 1._DP - - !- for OUTPUT hst/*.mtr.* - !%%% under benchmark %%% - domgdz = eps_r * ( sin(zz(iz)) & - + eps_hor * lmmq * sin( lmmq * zz(iz) - malpha ) & - + eps_mor * lmmqm1 * sin( lmmqm1 * zz(iz) - malpha ) & - + eps_por * lmmqp1 * sin( lmmqp1 * zz(iz) - malpha ) ) - domgdy = - eps_rnew / r_major * ( & - - ( sin( zz(iz) ) & - + eps_hor * lprd * sin( lmmq * zz(iz) - malpha ) & - + eps_mor * lprdm1 * sin( lmmqm1 * zz(iz) - malpha ) & - + eps_por * lprdp1 * sin( lmmqp1 * zz(iz) - malpha ) & - ) - (-1._DP/eps_r) * domgdz ) - domgdx = eps_rnew / r_major * ( & - - ( & - rdeps00 & - + rdeps1_0 * cos( zz(iz) ) & - + rdeps2_10 * cos( lmmq * zz(iz) - malpha ) & - + rdeps1_10 * cos( lmmqm1 * zz(iz) - malpha ) & - + rdeps3_10 * cos( lmmqp1 * zz(iz) - malpha ) & - + s_hat * zz(iz) * ( sin( zz(iz) ) & - + eps_hor * lprd * sin( lmmq * zz(iz) - malpha ) & - + eps_mor * lprdm1 * sin( lmmqm1 * zz(iz) - malpha ) & - + eps_por * lprdp1 * sin( lmmqp1 * zz(iz) - malpha ) ) & - ) - (-s_hat*zz(iz)/eps_r) * domgdz ) - gg(1,1) = 1._DP - gg(1,2) = s_hat*zz(iz) - gg(1,3) = 0._DP - gg(2,1) = gg(1,2) - gg(2,2) = 1._DP + (s_hat*zz(iz))**2 - gg(2,3) = 1._DP/(r_major*eps_r) - gg(3,1) = gg(1,3) - gg(3,2) = gg(2,3) - gg(3,3) = 1._DP/((r_major*eps_r)**2) - metric_l( 1,iz) = zz(iz) ! [ 1] - metric_l( 2,iz) = theta ! [ 2] - metric_l( 3,iz) = omg(iz) ! [ 3] - metric_l( 4,iz) = domgdx ! [ 4] - metric_l( 5,iz) = domgdy ! [ 5] - metric_l( 6,iz) = domgdz ! [ 6] - metric_l( 7,iz) = gg(1,1) ! [ 7] - metric_l( 8,iz) = gg(1,2) ! [ 8] - metric_l( 9,iz) = gg(1,3) ! [ 9] - metric_l(10,iz) = gg(2,2) ! [10] - metric_l(11,iz) = gg(2,3) ! [11] - metric_l(12,iz) = gg(3,3) ! [12] - metric_l(13,iz) = rootg(iz)! [13] - !------------------------- - -!!! for s-alpha !!! <--- the current version is the same as "analytic" - else if( trim(equib_type) == "s-alpha" .or. trim(equib_type) == "s-alpha-shift" ) then - - q_bar = q_0 - r_major = 1._DP ! in the R0 unit - - if (trim(equib_type) == "s-alpha") then - !--- s-alpha model without Shafranov shift - - alpha_MHD = 0._DP - else if (trim(equib_type) == "s-alpha-shift") then - !--- s-alpha model with Shafranov shift ---- - p_total = 0._DP - dp_totaldx = 0._DP - beta_total = 0._DP - do is = 0, ns-1 - p_total = p_total + fcs(is) * tau(is) / Znum(is) - dp_totaldx = dp_totaldx - fcs(is) * tau(is) / Znum(is) * (R0_Ln(is) + R0_Lt(is)) - beta_total = beta_total + 2._DP * beta * fcs(is) * tau(is) / Znum(is) - end do - alpha_MHD = - q_0**2 * r_major * beta_total * dp_totaldx / p_total - end if - - theta = zz(iz) - - omg(iz) = 1._DP - eps_r * cos( theta ) ! s-alpha with eps-expansion - !omg(iz) = 1._DP / (1._DP + eps_r * cos( theta )) ! for benchmark - - rootg(iz) = q_0*r_major/omg(iz) - dpara(iz) = dz* q_0 * r_major - - domgdz = eps_r * sin( theta ) - !domgdz = eps_r * sin( theta ) * omg(iz)**2 ! for benchmark - domgdx = -cos( theta ) / r_major - domgdy = 0._DP - - - gg(1,1) = 1._DP - gg(1,2) = s_hat*zz(iz) - alpha_MHD*sin(zz(iz)) ! with Shafranov shift - gg(1,3) = 0._DP - gg(2,1) = gg(1,2) - gg(2,2) = 1._DP + (s_hat*zz(iz) - alpha_MHD*sin(zz(iz)))**2 ! with Shafranov shift - gg(2,3) = 1._DP/(r_major*eps_r) - gg(3,1) = gg(1,3) - gg(3,2) = gg(2,3) - gg(3,3) = 1._DP/((r_major*eps_r)**2) - - kkx = -r_major * (q_0/q_bar) & - * ( gg(1,1)*gg(2,3) - gg(1,2)*gg(1,3) )*domgdz - kky = r_major * (q_bar/q_0) & - * ( domgdx - ( gg(1,2)*gg(2,3) - gg(2,2)*gg(1,3) )*domgdz ) - - do im = 0, nm - - vp(iz,im) = sqrt( 2._DP * mu(im) * omg(iz) ) - - mir(iz,im) = mu(im) * (q_0/q_bar) * domgdz / ( omg(iz)*rootg(iz) ) - - do iv = 1, 2*nv - !do my = ist_y, iend_y - ! do mx = -nx, nx - - !!!!kvd(mx,my,iz,iv,im) = & - !!!! ( vl(iv)**2 + omg(iz)*mu(im) ) / r_major & - !!!! * ( kkx*kx(mx) + kky*ky(my) ) & - !!!! * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) - - !!!!kvs(mx,my,iz,iv,im) = & - !!!! - sgn(ranks) * tau(ranks) / Znum(ranks) * ky(my) & - !!!! * ( R0_Ln(ranks) + R0_Lt(ranks) * ( 0.5_DP*vl(iv)**2 & - !!!! + omg(iz)*mu(im) - 1.5_DP ) ) & - !!!! * (q_bar/q_0) - - !%%% kvd = kx(mx) * vdx(iz,iv,im) + ky(my) * vdy(iz,iv,im) %%% - !%%% kvs = ky(my) * vsy(iz,iv,im) %%% - vdx(iz,iv,im) = & - ( vl(iv)**2 + omg(iz)*mu(im) ) / r_major & - * kkx & - * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) - - vdy(iz,iv,im) = & - ( vl(iv)**2 + omg(iz)*mu(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 ) ) & - * (q_bar/q_0) - - ! end do - !end do - end do - - end do ! im loop ends - - - ksq(:,:,iz) = 0._DP - do my = ist_y, iend_y - do mx = -nx, nx - ksq(mx,my,iz) = ( kx(mx) + ( s_hat * zz(iz) - alpha_MHD*sin(zz(iz)) ) & - * ky(my) )**2 + ky(my)**2 ! with Shafranov shift - end do - end do - - baxfactor = 1._DP - - !- for OUTPUT hst/*.mtr.* - - metric_l( 1,iz) = zz(iz) ! [ 1] - metric_l( 2,iz) = theta ! [ 2] - metric_l( 3,iz) = omg(iz) ! [ 3] - metric_l( 4,iz) = domgdx ! [ 4] - metric_l( 5,iz) = domgdy ! [ 5] - metric_l( 6,iz) = domgdz ! [ 6] - metric_l( 7,iz) = gg(1,1) ! [ 7] - metric_l( 8,iz) = gg(1,2) ! [ 8] - metric_l( 9,iz) = gg(1,3) ! [ 9] - metric_l(10,iz) = gg(2,2) ! [10] - metric_l(11,iz) = gg(2,3) ! [11] - metric_l(12,iz) = gg(3,3) ! [12] - metric_l(13,iz) = rootg(iz)! [13] - !------------------------- - - -!!! for circular MHD equilibrium !!! - else if( trim(equib_type) == "circ-MHD" ) then - - q_bar = dsqrt( 1._DP - eps_r**2 )*q_0 - r_major = 1._DP ! in the R0 unit - - theta = 2._DP*atan( sqrt( (1._DP+eps_r)/(1._DP-eps_r) ) & - * tan(zz(iz)/2._DP) ) - - omg(iz) = sqrt( q_bar**2 + eps_r**2 ) & - / ( 1._DP + eps_r*cos( theta ) ) / q_bar - - rootg(iz) = q_0*r_major*( 1._DP+eps_r*cos(theta) )**2 - dpara(iz) = dz * omg(iz) * rootg(iz) - - - domgdz = eps_r * sin(theta) * sqrt( q_bar**2 + eps_r**2 ) & - / ( 1._DP + eps_r * cos( theta ) )**2 & - / ( 1._DP - eps_r * cos( zz(iz)) ) / q_0 - - domgdx = -( cos(theta) & - - eps_r*(1._DP-s_hat+eps_r**2*q_0**2/q_bar**2) & - *(1._DP+eps_r*cos(theta))/(q_bar**2+eps_r**2) & - - eps_r*sin(theta)**2/(1._DP-eps_r**2) & - ) / ((1._DP + eps_r*cos(theta))**2) & - * sqrt(q_bar**2+eps_r**2) / q_bar / r_major - - domgdy = 0._DP - - gg(1,1) = (q_0/q_bar)**2 - gg(1,2) = ( s_hat*zz(iz)*q_0/q_bar - eps_r*sin(zz(iz))/(1._DP-eps_r**2) )*q_0/q_bar - gg(1,3) = - sin(zz(iz))/(1._DP-eps_r**2)/r_major*q_0/q_bar - gg(2,1) = gg(1,2) - gg(2,2) = (s_hat*zz(iz)*q_0/q_bar)**2 - 2._DP*q_0/q_bar*s_hat*zz(iz)*eps_r*sin(zz(iz))/(1._DP-eps_r**2) & - + (q_bar**2+eps_r**2)/((1._DP+eps_r*cos(theta))**2)/(q_0**2) & - + (eps_r*sin(zz(iz)))**2/(1._DP-eps_r**2)**2 - gg(2,3) = ( -s_hat*zz(iz)*q_0/q_bar*sin(zz(iz))/(1._DP-eps_r**2) & - + ((q_bar/q_0)**2)/((1._DP+eps_r*cos(theta))**2)/eps_r & - + eps_r*(sin(zz(iz))**2)/((1._DP-eps_r**2)**2) & - ) / r_major - gg(3,1) = gg(1,3) - gg(3,2) = gg(2,3) - gg(3,3) = ( ((q_bar/q_0)**2)/((1._DP+eps_r*cos(theta))**2)/(eps_r**2) & - + (sin(zz(iz))**2)/((1._DP-eps_r**2)**2) & - ) / (r_major**2) - - kkx = r_major*( -domgdy + (gg(1,3)*gg(1,2) - gg(1,1)*gg(2,3))*domgdz/omg(iz)**2 ) - kky = r_major*( domgdx + (gg(1,3)*gg(2,2) - gg(1,2)*gg(2,3))*domgdz/omg(iz)**2 ) - - do im = 0, nm - - vp(iz,im) = sqrt( 2._DP * mu(im) * omg(iz) ) - - mir(iz,im) = mu(im) * domgdz / ( omg(iz)*rootg(iz) ) - - do iv = 1, 2*nv - !do my = ist_y, iend_y - ! do mx = -nx, nx - - !!!!kvd(mx,my,iz,iv,im)= & - !!!! ( vl(iv)**2 + omg(iz)*mu(im) ) / ( r_major*omg(iz) ) & - !!!! * ( kkx*kx(mx) + kky*ky(my) ) & - !!!! * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) - - !!!!kvs(mx,my,iz,iv,im) = & - !!!! - sgn(ranks) * tau(ranks) / Znum(ranks) * ky(my) & - !!!! * ( R0_Ln(ranks) + R0_Lt(ranks) * ( 0.5_DP*vl(iv)**2 & - !!!! + omg(iz)*mu(im) - 1.5_DP ) ) - - !%%% kvd = kx(mx) * vdx(iz,iv,im) + ky(my) * vdy(iz,iv,im) %%% - !%%% kvs = ky(my) * vsy(iz,iv,im) %%% - vdx(iz,iv,im)= & - ( vl(iv)**2 + omg(iz)*mu(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) ) & - * 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 ) ) - - - ! end do - !end do - end do - - end do ! im loop ends - - - ksq(:,:,iz) = 0._DP - do my = ist_y, iend_y - do mx = -nx, nx - ksq(mx,my,iz) = (kx(mx)**2)*gg(1,1) & - + 2._DP*kx(mx)*ky(my)*gg(1,2) & - + (ky(my)**2)*gg(2,2) - end do - end do - - baxfactor = 1._DP - - !- for OUTPUT hst/*.mtr.* - - metric_l( 1,iz) = zz(iz) ! [ 1] - metric_l( 2,iz) = theta ! [ 2] - metric_l( 3,iz) = omg(iz) ! [ 3] - metric_l( 4,iz) = domgdx ! [ 4] - metric_l( 5,iz) = domgdy ! [ 5] - metric_l( 6,iz) = domgdz ! [ 6] - metric_l( 7,iz) = gg(1,1) ! [ 7] - metric_l( 8,iz) = gg(1,2) ! [ 8] - metric_l( 9,iz) = gg(1,3) ! [ 9] - metric_l(10,iz) = gg(2,2) ! [10] - metric_l(11,iz) = gg(2,3) ! [11] - metric_l(12,iz) = gg(3,3) ! [12] - metric_l(13,iz) = rootg(iz)! [13] - !------------------------- - -!!!! for VMEC equilibrium !!! -! else if( trim(equib_type) == "vmec" ) then -! -! q_bar = q_0 -! theta = zz(iz) -! -! call vmecin_coeff( rad_a, R0_unit, rho2R_0, q_input, theta, & -! alpha_fix, r_0, r_minor, s_hat, & -! gdwss, gdwtt, gdwzz, gdwst, gdwsz, gdwtz, & -! gupss, guptt, gupzz, gupst, gupsz, guptz, & -! babs, Bs, Bth, Bzt, dBds, dBdt, dBdz, & -! dBdt_mir, vmec_rootg, rootgft, rootgbz ) -! -! omg(iz) = babs -! -! rootg(iz) = vmec_rootg * R0_unit * R0_unit * R0_unit -! dpara(iz) = dz * babs * rootgft * b0b00 -! -! -! -! do im = 0, nm -! -! vp(iz,im) = sqrt( 2._DP * mu(im) * omg(iz) ) -! mir(iz,im) = mu(im) * dBdt_mir / babs / rootgft / b0b00 -! -! do iv = 1, 2*nv -! do my = ist_y, iend_y -! do mx = -nx, nx -! -! kvd(mx,my,iz,iv,im) = & -! - (( vl(iv)**2 + omg(iz)*mu(im) ) / rootgbz /babs/babs/babs ) & -! * ((r_0/q_0) * ky(my) & -! * ( ( (Bs/r_a) + Bzt * (q_0/r_0) * s_hat * zz(iz) ) * dBdt & -! +( (Bs/r_a) * q_0 - Bth * (q_0/r_0) * s_hat * zz(iz) ) * dBdz & -! -( Bth + Bzt * q_0 ) * dBds / r_a ) & -! + kx(mx) * ( Bzt * dBdt - Bth * dBdz )) & -! * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) -! ! --- k*v_d term -! -! kvs(mx,my,iz,iv,im) = & -! - sgn(ranks) * ky(my) & -! * ((r_0/q_0) * (Bth + Bzt * q_0) / rootgbz / babs / babs) & -! * ( R0_Ln(ranks) & -! + R0_Lt(ranks) * (0.5_DP*vl(iv)**2 + omg(iz)*mu(im) - 1.5_DP) ) & -! * tau(ranks) / Znum(ranks) -! ! --- k*v_* term -! end do -! end do -! end do -! -! end do ! im loop ends -! -! -! do my = ist_y, iend_y -! do mx = -nx, nx -! ksq(mx,my,iz) = (r_a * kx(mx))**2 * gupss & -! + ky(my)**2 * ( (r_0/q_0)**2 & -! * ( gupzz + guptt * q_0 **2 - guptz * 2._DP * q_0 ) & -! + 2._DP * s_hat * (r_0/q_0) * zz(iz) * r_a * ( gupst * q_0 - gupsz ) & -! + r_a * r_a * gupss * (s_hat**2) * (zz(iz)**2) ) & -! + (r_a * kx(mx)) * ky(my) * 2._DP * ( (r_0/q_0) & -! * ( gupst * q_0 - gupsz ) + r_a * gupss * s_hat * zz(iz) ) -! ! --- squere of k_perp -! end do -! end do -! -! baxfactor = b0b00 ! --- For the use in caldlt -! - -! this is new vmec-BoozXform interface by M. Nakata & M. Nunami (Aug. 2016) - else if( trim(equib_type) == "vmec" ) then - - q_bar = q_0 - isw = 1 - r_major = 1._DP ! in the R0 unit - - call vmecbzx_boozx_coeff( isw, nss, ntheta, nzeta, s_input, iz, zz(iz), lz_l, & ! input - s_0, q_0, s_hat, eps_r, phi_ax, & ! output - omg(iz), rootg(iz), domgdx, domgdz, domgdy, & - gg(1,1), gg(1,2), gg(1,3), gg(2,2), & - gg(2,3), gg(3,3) ) - - dpara(iz) = dz * omg(iz) * rootg(iz) - - kkx = r_major*( -domgdy + (gg(1,3)*gg(1,2) - gg(1,1)*gg(2,3))*domgdz/omg(iz)**2 ) - kky = r_major*( domgdx + (gg(1,3)*gg(2,2) - gg(1,2)*gg(2,3))*domgdz/omg(iz)**2 ) - - do im = 0, nm - - vp(iz,im) = sqrt( 2._DP * mu(im) * omg(iz) ) - - mir(iz,im) = mu(im) * domgdz / ( omg(iz)*rootg(iz) ) - - do iv = 1, 2*nv - !do my = ist_y, iend_y - ! do mx = -nx, nx - - !!!!kvd(mx,my,iz,iv,im) = & - !!!! ( vl(iv)**2 + omg(iz)*mu(im) ) / ( r_major*omg(iz) ) & - !!!! * ( kkx*kx(mx) + kky*ky(my) ) & - !!!! * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) & - !!!! - real(ibprime,kind=DP) * vl(iv)**2 / r_major / omg(iz)**2 & ! grad-p (beta-prime) term - !!!! * ( beta*(R0_Ln(ranks) + R0_Lt(ranks))*ky(my) ) & - !!!! * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) - - !!!!kvs(mx,my,iz,iv,im) = & - !!!! - sgn(ranks) * tau(ranks) / Znum(ranks) * ky(my) & - !!!! * ( R0_Ln(ranks) + R0_Lt(ranks) * ( 0.5_DP*vl(iv)**2 & - !!!! + omg(iz)*mu(im) - 1.5_DP ) ) - - !%%% kvd = kx(mx) * vdx(iz,iv,im) + ky(my) * vdy(iz,iv,im) %%% - !%%% kvs = ky(my) * vsy(iz,iv,im) %%% - vdx(iz,iv,im) = & - ( vl(iv)**2 + omg(iz)*mu(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) ) & - * kky & - * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) & - - real(ibprime,kind=DP) * vl(iv)**2 / r_major / omg(iz)**2 & ! grad-p (beta-prime) term - * ( beta*(R0_Ln(ranks) + R0_Lt(ranks)) ) & - * ( 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 ) ) - - - ! end do - !end do - end do - - end do ! im loop ends - - - ksq(:,:,iz) = 0._DP - do my = ist_y, iend_y - do mx = -nx, nx - ksq(mx,my,iz) = (kx(mx)**2)*gg(1,1) & - + 2._DP*kx(mx)*ky(my)*gg(1,2) & - + (ky(my)**2)*gg(2,2) - end do - end do - - baxfactor = 1._DP - - !- for OUTPUT hst/*.mtr.* - - metric_l( 1,iz) = zz(iz) ! [ 1] - metric_l( 2,iz) = phi_ax ! [ 2] Axisymetric toroidal angle - metric_l( 3,iz) = omg(iz) ! [ 3] - metric_l( 4,iz) = domgdx ! [ 4] - metric_l( 5,iz) = domgdy ! [ 5] - metric_l( 6,iz) = domgdz ! [ 6] - metric_l( 7,iz) = gg(1,1) ! [ 7] - metric_l( 8,iz) = gg(1,2) ! [ 8] - metric_l( 9,iz) = gg(1,3) ! [ 9] - metric_l(10,iz) = gg(2,2) ! [10] - metric_l(11,iz) = gg(2,3) ! [11] - metric_l(12,iz) = gg(3,3) ! [12] - metric_l(13,iz) = rootg(iz)! [13] - !------------------------- - - - else if( trim(equib_type) == "eqdsk" ) then - - q_bar = q_0 - isw = 1 - r_major = 1._DP ! in the R0 unit - - call igs_coeff( isw, mc_type, nss, ntheta, s_input, zz(iz), lz_l, & ! input - s_0, q_0, s_hat, eps_r, theta, & ! output - omg(iz), rootg(iz), domgdx, domgdz, domgdy, & - gg(1,1), gg(1,2), gg(1,3), gg(2,2), & - gg(2,3), gg(3,3) ) - - dpara(iz) = dz * omg(iz) * rootg(iz) - - kkx = r_major*( -domgdy + (gg(1,3)*gg(1,2) - gg(1,1)*gg(2,3))*domgdz/omg(iz)**2 ) - kky = r_major*( domgdx + (gg(1,3)*gg(2,2) - gg(1,2)*gg(2,3))*domgdz/omg(iz)**2 ) - - do im = 0, nm - - vp(iz,im) = sqrt( 2._DP * mu(im) * omg(iz) ) - - mir(iz,im) = mu(im) * domgdz / ( omg(iz)*rootg(iz) ) - - do iv = 1, 2*nv - !do my = ist_y, iend_y - ! do mx = -nx, nx - - !!!!kvd(mx,my,iz,iv,im) = & - !!!! ( vl(iv)**2 + omg(iz)*mu(im) ) / ( r_major*omg(iz) ) & - !!!! * ( kkx*kx(mx) + kky*ky(my) ) & - !!!! * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) & - !!!! - real(ibprime,kind=DP) * vl(iv)**2 / r_major / omg(iz)**2 & ! grad-p (beta-prime) term - !!!! * ( beta*(R0_Ln(ranks) + R0_Lt(ranks))*ky(my) ) & - !!!! * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) - - !!!!kvs(mx,my,iz,iv,im) = & - !!!! - sgn(ranks) * tau(ranks) / Znum(ranks) * ky(my) & - !!!! * ( R0_Ln(ranks) + R0_Lt(ranks) * ( 0.5_DP*vl(iv)**2 & - !!!! + omg(iz)*mu(im) - 1.5_DP ) ) - - !%%% kvd = kx(mx) * vdx(iz,iv,im) + ky(my) * vdy(iz,iv,im) %%% - !%%% kvs = ky(my) * vsy(iz,iv,im) %%% - vdx(iz,iv,im) = & - ( vl(iv)**2 + omg(iz)*mu(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) ) & - * kky & - * ( sgn(ranks) * tau(ranks) / Znum(ranks) ) & - - real(ibprime,kind=DP) * vl(iv)**2 / r_major / omg(iz)**2 & ! grad-p (beta-prime) term - * ( beta*(R0_Ln(ranks) + R0_Lt(ranks)) ) & - * ( 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 ) ) - - - ! end do - !end do - end do - - end do ! im loop ends - - - ksq(:,:,iz) = 0._DP - do my = ist_y, iend_y - do mx = -nx, nx - ksq(mx,my,iz) = (kx(mx)**2)*gg(1,1) & - + 2._DP*kx(mx)*ky(my)*gg(1,2) & - + (ky(my)**2)*gg(2,2) - end do - end do - - baxfactor = 1._DP - - !- for OUTPUT hst/*.mtr.* - - metric_l( 1,iz) = zz(iz) ! [ 1] - metric_l( 2,iz) = theta ! [ 2] - metric_l( 3,iz) = omg(iz) ! [ 3] - metric_l( 4,iz) = domgdx ! [ 4] - metric_l( 5,iz) = domgdy ! [ 5] - metric_l( 6,iz) = domgdz ! [ 6] - metric_l( 7,iz) = gg(1,1) ! [ 7] - metric_l( 8,iz) = gg(1,2) ! [ 8] - metric_l( 9,iz) = gg(1,3) ! [ 9] - metric_l(10,iz) = gg(2,2) ! [10] - metric_l(11,iz) = gg(2,3) ! [11] - metric_l(12,iz) = gg(3,3) ! [12] - metric_l(13,iz) = rootg(iz)! [13] - !------------------------- - - - else - - write( olog, * ) " # wrong choice of the equilibrium " - call flush(olog) - call MPI_Finalize(ierr_mpi) - stop - - end if - - - 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) ) & - * 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) ) - call math_j2( kmo, j2(mx,my,iz,im) ) - end do - end do - end do - - - do my = ist_y, iend_y - do mx = -nx, nx - bb = ksq(mx,my,iz) / omg(iz)**2 & - * tau(ranks)*Anum(ranks)/(Znum(ranks)**2) - call math_g0( bb, g0(mx,my,iz) ) - end do - end do - - -!!!! debug (Jan 2012) -! write( olog, fmt="(1p,10e15.7)" ) & -! zz(iz), omg(iz), mir(iz,0), dpara(iz), jcob(iz), & -! ksq(1,2,iz), kvs(1,2,iz,1,0), kvd(1,2,iz,1,0), j0(1,2,iz,0), g0(1,2,iz) -!!!! debug (Jan 2012) - - - end do ! iz loop ends - -!- OUTPUT ascii data hst/*.mtr.* - - call MPI_gather(metric_l(1,-nz), num_omtr*2*nz, MPI_DOUBLE_PRECISION, & - metric_g(1,-global_nz), num_omtr*2*nz, MPI_DOUBLE_PRECISION, & - 0, zsp_comm_world, ierr_mpi) - if ( rankg == 0 ) then - do iz = -global_nz, global_nz-1 - write( omtr, fmt="(f15.8,SP,256E24.14e3)") metric_g(:,iz) - end do - call flush(omtr) - end if -!--------------------------------- - -! --- operator settings --- - - - cfsrf = 0._DP - cfsrf_l = 0._DP - do iz = -nz, nz-1 -! cfsrf_l = cfsrf_l + 1._DP / omg(iz) - cfsrf_l = cfsrf_l + rootg(iz) - ! normalization coefficient for - ! the surface average - end do - - call MPI_Allreduce( cfsrf_l, cfsrf, 1, MPI_DOUBLE_PRECISION, & - MPI_SUM, zsp_comm_world, ierr_mpi ) - - - ! --- debug - ! write( olog, * ) " *** z, omg " - ! do iz = -nz, nz-1 - ! write( olog, * ) zz(iz), omg(iz) - ! end do - ! write( olog, * ) - - - if ( vel_rank == 0 ) then - do iz = -nz, nz-1 - !dvp(iz) = vp(iz,1) - dvp(iz) = sqrt( 2._DP * (0.5_DP * dm**2) * omg(iz) ) - end do - end if - - call MPI_Bcast( dvp, 2*nz, MPI_DOUBLE_PRECISION, 0, & - vel_comm_world, ierr_mpi ) - - - do my = ist_y_g, iend_y_g - ck(my-ist_y_g) = exp( ui * 2._DP * pi * del_c & - * real( n_tht * my, kind=DP ) ) - dj(my-ist_y_g) = - m_j * n_tht * my - ! del_c = q_0*n_alp-int(q_0*n_alp) - ! m_j = 2*n_alp*q_d - end do - - - 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) ) & - / sqrt( twopi**3 ) - end do - end do - end do - - allocate( ww(-nx:nx,0:ny,-nz:nz-1) ) - -! --- GK polarization factor for efield calculation - fct_poisson(:,:,:) = 0._DP - fct_e_energy(:,:,:) = 0._DP - - ww(:,:,:) = 0._DP - do iz = -nz, nz-1 - do my = ist_y, iend_y - do mx = -nx, nx - - if ( rankw == 0 .and. mx == 0 .and. my == 0 ) then !- (0,0) mode - - fct_poisson(mx,my,iz) = 0._DP - fct_e_energy(mx,my,iz) = 0._DP - - else - - ww(mx,my,iz) = lambda_i * ksq(mx,my,iz) - do is = 0, ns-1 - bb = ksq(mx,my,iz) / omg(iz)**2 & - * tau(is)*Anum(is)/(Znum(is)**2) - call math_g0( bb, gg0 ) - ww(mx,my,iz) = ww(mx,my,iz) & - + Znum(is) * fcs(is) / tau(is) * ( 1._DP - gg0 ) - end do - fct_poisson(mx,my,iz) = 1._DP / ww(mx,my,iz) - fct_e_energy(mx,my,iz) = ww(mx,my,iz) - - end if - - end do - end do - end do - - -! --- ZF-factor for adiabatic model - if ( ns == 1 ) then - - ww(:,:,:) = 0._DP - do iz = -nz, nz-1 - my = 0 - do mx = -nx, nx - ww(mx,my,iz) = ( 1._DP - g0(mx,my,iz) ) & - / ( 1._DP - g0(mx,my,iz) + tau(0)*tau_ad ) - end do - end do - - call intgrl_fsrf ( ww, fctgt ) - - if ( rankw == 0 ) then - fctgt(0) = ( 1._DP - g0(0,0,0) ) / ( 1._DP - g0(0,0,0) + tau(0)*tau_ad ) - ! g0(0,0,iz) has no z dependence - endif - - endif - - deallocate( ww ) - - allocate( wf(-nx:nx,0:ny,-nz:nz-1,1:2*nv,0:nm) ) - allocate( nw(-nx:nx,0:ny,-nz:nz-1) ) - wf(:,:,:,:,:) = ( 0._DP, 0._DP ) - nw(:,:,:) = ( 0._DP, 0._DP ) - -! --- GK polarization factor for mfield calculation - fct_ampere(:,:,:) = 0._DP - fct_m_energy(:,:,:) = 0._DP - - if ( beta .ne. 0._DP ) 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 - wf(mx,my,iz,iv,im) = Znum(ranks) * fcs(ranks) / Anum(ranks) & - * vl(iv)**2 * j0(mx,my,iz,im)**2 * fmx(iz,iv,im) - end do - end do - end do - end do - end do - - call intgrl_v0_moment_ms ( wf, nw ) - - do iz = -nz, nz-1 - do my = ist_y, iend_y - do mx = -nx, nx - fct_ampere(mx,my,iz) = 1._DP / real( ksq(mx,my,iz) + beta * nw(mx,my,iz), kind=DP ) - fct_m_energy(mx,my,iz) = ksq(mx,my,iz) / beta - end do - end do - end do - - if ( rankw == 0 ) then - do iz = -nz, nz-1 - fct_ampere(0,0,iz) = 0._DP - fct_m_energy(0,0,iz) = 0._DP - end do - end if - - end if - - deallocate( wf ) - deallocate( nw ) +! --- read GKV namelist relating to configurations --- + call geom_read_nml +! --- coordinate settings (time-indep.) --- + call geom_init_kxkyzvm(lx, ly, eps_r) ! --- set collision frequencies and v-space functions for multi-species GK collision read(inml,nml=nu_ref) call colli_set_param(q_0, eps_r, nust) - !if (trim(time_advnc) == "imp_colli" .or. trim(time_advnc) == "auto_init") then - if (trim(col_type) == "full" .or. trim(col_type) == "lorentz" .or. trim(time_advnc) == "imp_colli") then - call colliimp_set_param - end if - !!! call colliimp_set_param - call dtc_init( lx, ly, vmax ) - write( olog, * ) " # Collision parameters" write( olog, * ) "" write( olog, * ) " # Nref [m^-3] = ", Nref @@ -1826,6 +444,17 @@ SUBROUTINE set_cnfig write( olog, * ) +! --- coordinate settings (explicitly time-dependent metrics) --- + call geom_init_metric + +! --- operator settings (time-dependent through metrics) --- + call geom_set_operators + if (trim(col_type) == "full" .or. trim(col_type) == "lorentz" .or. trim(time_advnc) == "imp_colli") then + call colliimp_set_param + end if + +! --- initial estimate of time steps --- + call dtc_init( lx, ly, vmax ) END SUBROUTINE set_cnfig @@ -1969,6 +598,15 @@ SUBROUTINE set_value( ff, phi, Al, hh, time ) end if + !%%% For shearflow rotating flux tube model %%% + if (gamma_e /= 0._DP .and. trim(flag_shearflow) =="rotating") then + call geom_reset_time(time) + if (trim(col_type) == "full" .or. trim(col_type) == "lorentz" .or. trim(time_advnc) == "imp_colli") then + call colliimp_set_param + end if + end if + !%%% + call bndry_zvm_bound_f( ff ) call fld_esfield( ff, phi ) diff --git a/src/gkvp_vmecbzx.f90 b/src/gkvp_vmecbzx.f90 index f2c9f98..25bbb64 100644 --- a/src/gkvp_vmecbzx.f90 +++ b/src/gkvp_vmecbzx.f90 @@ -6,6 +6,9 @@ MODULE GKV_vmecbzx ! ! Update history of gkvp_vmecbxz.f90 ! -------------- +! gkvp_f0.62 (S. Maeyama, Mar 2023) +! - Input of vmecbzx_boozxcoef is modified from local iz for each rankz +! to global index giz. ! gkvp_f0.57 (S. Maeyama, Oct 2020) ! - Version number f0.57 is removed from filename. ! @@ -43,9 +46,9 @@ SUBROUTINE vmecbzx_boozx_read( nss, ntheta, nzeta ) implicit none integer, intent(in) :: nss, ntheta, nzeta - integer :: is, jj, ierr, ibzx + integer :: ibzx character(512) :: f_bozx - character(512) :: env_string !fj +! character(512) :: env_string !fj namelist /bozxf/ f_bozx @@ -132,14 +135,20 @@ END SUBROUTINE vmecbzx_boozx_read !---------------------------------------------------------------------------------- - SUBROUTINE vmecbzx_boozx_coeff( isw, nss, ntheta, nzeta, s_input, iz, zz, lz_l, & ! input +!smae start 202303 +! SUBROUTINE vmecbzx_boozx_coeff( isw, nss, ntheta, nzeta, s_input, iz, zz, lz_l, & ! input + SUBROUTINE vmecbzx_boozx_coeff( isw, nss, ntheta, nzeta, s_input, giz, zz, lz_l, & ! input +!smae end 202303 s_0, q_0, s_hat, eps_r, phi_ax, & ! output omg, rootg, domgdx, domgdz, domgdy, & gg11, gg12, gg13, gg22, & gg23, gg33 ) !---------------------------------------------------------------------------------- - integer, intent(in) :: isw, nss, ntheta, nzeta, iz +!smae start 202303 +! integer, intent(in) :: isw, nss, ntheta, nzeta, iz + integer, intent(in) :: isw, nss, ntheta, nzeta, giz +!smae end 202303 real(kind=DP), intent(in) :: s_input, zz, lz_l real(kind=DP), intent(inout) :: s_0, q_0, s_hat, eps_r, phi_ax @@ -147,8 +156,8 @@ SUBROUTINE vmecbzx_boozx_coeff( isw, nss, ntheta, nzeta, s_input, iz, zz, real(kind=DP), intent(out) :: gg11, gg12, gg13, gg22, gg23, gg33 ! --- local variables - integer :: is, jj, is0, iz0, nz0, jj0, zt0, giz - real(kind=DP) :: zz0, eps_a + integer :: is0, jj0, zt0 + real(kind=DP) :: eps_a ! is0 = nint(s_input*(nss+1)) @@ -164,7 +173,9 @@ SUBROUTINE vmecbzx_boozx_coeff( isw, nss, ntheta, nzeta, s_input, iz, zz, else if ( isw == 1 ) then s_0 = ss_bz(is0) - giz = (-global_nz + 2*nz*rankz + iz + nz ) +!smae start 202303 +! giz = (-global_nz + 2*nz*rankz + iz + nz ) +!smae end 202303 jj0 = giz + ntheta/2