From becdad2c35720aa4938f780994789835fe29397d Mon Sep 17 00:00:00 2001 From: lazersos Date: Thu, 22 Aug 2024 14:41:30 +0200 Subject: [PATCH 01/47] PENTA: Added wrapper main function so PENTA can be called like a subroutine. --- PENTA/ObjectList | 1 + PENTA/PENTA.dep | 3 + PENTA/Sources/penta_imp.f | 143 +++++++++---------------------------- PENTA/Sources/penta_main.f | 126 ++++++++++++++++++++++++++++++++ 4 files changed, 163 insertions(+), 110 deletions(-) create mode 100644 PENTA/Sources/penta_main.f diff --git a/PENTA/ObjectList b/PENTA/ObjectList index 488576fbf..6a0b86914 100644 --- a/PENTA/ObjectList +++ b/PENTA/ObjectList @@ -1,4 +1,5 @@ ObjectFiles = \ +penta_main.o \ penta_functions_mod.o \ penta_modules.o \ quanc8.o \ diff --git a/PENTA/PENTA.dep b/PENTA/PENTA.dep index db129c506..960e6f136 100644 --- a/PENTA/PENTA.dep +++ b/PENTA/PENTA.dep @@ -1,4 +1,7 @@ # Microsoft Developer Studio Generated Dependency File, included by DKES.mak +penta_main.o : \ + penta_modules.o + penta_imp.o : \ read_input_file_mod.o \ diff --git a/PENTA/Sources/penta_imp.f b/PENTA/Sources/penta_imp.f index b90c0b0bc..b4acb8633 100644 --- a/PENTA/Sources/penta_imp.f +++ b/PENTA/Sources/penta_imp.f @@ -1,90 +1,5 @@ -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c -c -PENTA with Impurities- -c -ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c To run,type: PENTA_imp [command line arguments] -c -c where the command line arguments are: -c -c 1 unique identifier text on DKES database files -c i.e., for (m,l,n)star_lijs_***, this would be -c the text that replaces *** for the specific surface (ex: hsx_s10) -c 2 Er,min -c 3 Er,max - args 2 and 3 specify the interval in Er over which -c the search is done for the ambipolar electric field root. -c Er,min and Er,max are dimensionless normalized electric -c field variables, i.e.: eEr,min/max/kT_e -c - Note if the logical variable input_is_Er below is true then -c the search interval is set by args 2 and 3 as Er (V/cm) -c 4 flux surface number -c 5 parameter that determines whether output data should be -c written into new files (=0) or appended to -c existing files (=1). Typically,when a script is run for a -c sequence of flux surfaces, this is set to 0 for the first -c surface and to 1 for subsequent surfaces -c 6 extension of the profile_data_*** file (i.e., the *** text) -c that comes from extracting the appropriate profiles from -c the VMEC run using the profile_extractor.f code -c 7 unique identifier on plasma profiles file, i.e., the *** -c text in the filename plasma_profiles***.dat. This allows -c multiple profiles to be run without having to rename files -c 8 in T*V/m -c -c -c Required input files: -c -c [l,m,n]star_lijs_***_s# - where *** is arg1 above and # arg 4. -c These files contain the values of efield and cmul and the -c corresponding normalized monoenergetic coefficients. Note -c that L* should be L*/e**2 where e is the elementary charge. -c -c ion_params - A file containing the namelist "ion_params" which -c defines the variables num_ion_species, Z_ion_init and -c miomp_init which are the number of ion species (integer), -c the corresponding ion charge numbers (real, array) and -c the ion to proton mass ratio (real, array) for each -c non-electron species. -c -c profile_data_*** - where *** is arg1 above. Contains geometry -c variables, most likely from a VMEC run. -c -c plasma_profilesXXX.dat - A file containing the plasma profile -c data, where XXX is arg8 above. This file has the format: -c row 1: number of radial points for which data is specified. -c all other rows: r/a, ne, Te, Ti1, ni1, Ti2, ni2 ... -c where the densities are in units of 10**12/cc and the -c temperatures in eV and i1, i2... are the ion species. -c -c Utilde2_profile - Contains the quantity , where U is the -c PS flow function as defined in Sugama and Nishimura. The -c first row is the number of points, then r/a points and the -c corresponding value. -c -c Output files: -c -c "fluxes_vs_roa", "fluxes_vs_Er", "flows_vs_roa", "plasma_profiles_check" -c -c -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c NOTES -c - Be sure to check parameters below -c - based on Sugama, Nishimura PoP 9 4637 (2002) and -c Sugama, Nishimura PoP 15, 042502 (2008) -c - Some code used from orignal PENTA code, written by Don Spong and -c modified by Jeremy Lore -c -c 7/2009 JL -ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c PROGRAMMING NOTES -c -speed optimization -c -recursive ambipolar solver -c -fix numargs -ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - program penta_imp + SUBROUTINE penta_imp(coeff_ext, Er_min, Er_max, js, i_append, + 1 run_ident, pprof_char, B_Eprl_in) use penta_kind_mod use phys_const use read_input_file_mod @@ -107,12 +22,21 @@ program penta_imp integer(iknd),parameter :: numKsteps_set = 1000_iknd !Number of K points (linear) for convolution (used if use_quanc8_set=.false.) logical, parameter :: log_interp_set = .true. !If true, logarithmic interpolation of DKES coefficients is used - integer(iknd) :: num_species, iroot, js, i_append, + character(100), intent(inout) :: coeff_ext + real(rknd), intent(inout) :: Er_min + real(rknd), intent(inout) :: Er_max + integer(iknd), intent(inout) :: js + integer(iknd), intent(inout) :: i_append + character(100), intent(inout) :: run_ident + character(10), intent(inout) :: pprof_char + real(rknd), intent(inout), OPTIONAL :: B_Eprl_in + + integer(iknd) :: num_species, iroot, 1 num_ion_species, inv_err, nn_inv, j, ispec, tmp_ind, ie, 2 ispec1, ispec2, spec1_ind, spec2_ind, numargs real(rknd) :: vth_e, loglambda, sigma_S, Er_test, U2, - 1 lX_sum1, lX_sum2, J_E_tot, Er_min, Er_max, B_Eprl, J_E_cl, + 1 lX_sum1, lX_sum2, J_E_tot, J_E_cl, B_Eprl, 2 eaEr_o_kTe real(rknd), dimension(:), allocatable :: Z_ion,ion_mass,vth_i, @@ -130,9 +54,8 @@ program penta_imp real(rknd), dimension(num_ion_max) :: Z_ion_init,miomp_init character(1) :: tb = char(9) - character(10) :: pprof_char character(100) :: arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, - 1 arg9, coeff_ext, run_ident, format_tmp, str_num, fstatus, fpos + 1 arg9, format_tmp, str_num, fstatus, fpos namelist / ion_params / num_ion_species,Z_ion_init,miomp_init @@ -154,6 +77,7 @@ program penta_imp !default values B_Eprl=0._rknd + IF (PRESENT(B_Eprl_in)) B_Eprl= B_Eprl_in !Read namelist file for ion parameters open(iu_nl,file="ion_params",status="old",form="formatted") @@ -162,25 +86,25 @@ program penta_imp !Get command line arguments !numargs = iargc() - numargs = 8 - call getarg(1, arg1) - call getarg(2, arg2) - call getarg(3, arg3) - call getarg(4, arg4) - call getarg(5, arg5) - call getarg(6, arg6) - call getarg(7, arg7) - if (numargs .eq. 8) call getarg(8, arg8) +! numargs = 8 +! call getarg(1, arg1) +! call getarg(2, arg2) +! call getarg(3, arg3) +! call getarg(4, arg4) +! call getarg(5, arg5) +! call getarg(6, arg6) +! call getarg(7, arg7) +! if (numargs .eq. 8) call getarg(8, arg8) !store command line args - coeff_ext = trim(adjustl(arg1)) - read(arg2,*) Er_min - read(arg3,*) Er_max - read(arg4,*) js - read(arg5,*) i_append - run_ident = trim(adjustl(arg6)) - pprof_char = trim(adjustl(arg7)) - if (numargs .eq. 8) read(arg8,*) B_Eprl +! coeff_ext = trim(adjustl(arg1)) +! read(arg2,*) Er_min +! read(arg3,*) Er_max +! read(arg4,*) js +! read(arg5,*) i_append +! run_ident = trim(adjustl(arg6)) +! pprof_char = trim(adjustl(arg7)) +! if (numargs .eq. 8) read(arg8,*) B_Eprl !allocate variables by number of species defined num_species=num_ion_species+1_iknd @@ -509,5 +433,4 @@ program penta_imp close(iu_fvEr_out) close(iu_flows_out) - end program penta_imp - + END SUBROUTINE penta_imp \ No newline at end of file diff --git a/PENTA/Sources/penta_main.f b/PENTA/Sources/penta_main.f new file mode 100644 index 000000000..541f6a2dd --- /dev/null +++ b/PENTA/Sources/penta_main.f @@ -0,0 +1,126 @@ +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c +c -PENTA with Impurities- +c +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c To run,type: PENTA_imp [command line arguments] +c +c where the command line arguments are: +c +c 1 unique identifier text on DKES database files +c i.e., for (m,l,n)star_lijs_***, this would be +c the text that replaces *** for the specific surface (ex: hsx_s10) +c 2 Er,min +c 3 Er,max - args 2 and 3 specify the interval in Er over which +c the search is done for the ambipolar electric field root. +c Er,min and Er,max are dimensionless normalized electric +c field variables, i.e.: eEr,min/max/kT_e +c - Note if the logical variable input_is_Er below is true then +c the search interval is set by args 2 and 3 as Er (V/cm) +c 4 flux surface number +c 5 parameter that determines whether output data should be +c written into new files (=0) or appended to +c existing files (=1). Typically,when a script is run for a +c sequence of flux surfaces, this is set to 0 for the first +c surface and to 1 for subsequent surfaces +c 6 extension of the profile_data_*** file (i.e., the *** text) +c that comes from extracting the appropriate profiles from +c the VMEC run using the profile_extractor.f code +c 7 unique identifier on plasma profiles file, i.e., the *** +c text in the filename plasma_profiles***.dat. This allows +c multiple profiles to be run without having to rename files +c 8 in T*V/m +c +c +c Required input files: +c +c [l,m,n]star_lijs_***_s# - where *** is arg1 above and # arg 4. +c These files contain the values of efield and cmul and the +c corresponding normalized monoenergetic coefficients. Note +c that L* should be L*/e**2 where e is the elementary charge. +c +c ion_params - A file containing the namelist "ion_params" which +c defines the variables num_ion_species, Z_ion_init and +c miomp_init which are the number of ion species (integer), +c the corresponding ion charge numbers (real, array) and +c the ion to proton mass ratio (real, array) for each +c non-electron species. +c +c profile_data_*** - where *** is arg1 above. Contains geometry +c variables, most likely from a VMEC run. +c +c plasma_profilesXXX.dat - A file containing the plasma profile +c data, where XXX is arg8 above. This file has the format: +c row 1: number of radial points for which data is specified. +c all other rows: r/a, ne, Te, Ti1, ni1, Ti2, ni2 ... +c where the densities are in units of 10**12/cc and the +c temperatures in eV and i1, i2... are the ion species. +c +c Utilde2_profile - Contains the quantity , where U is the +c PS flow function as defined in Sugama and Nishimura. The +c first row is the number of points, then r/a points and the +c corresponding value. +c +c Output files: +c +c "fluxes_vs_roa", "fluxes_vs_Er", "flows_vs_roa", "plasma_profiles_check" +c +c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c NOTES +c - Be sure to check parameters below +c - based on Sugama, Nishimura PoP 9 4637 (2002) and +c Sugama, Nishimura PoP 15, 042502 (2008) +c - Some code used from orignal PENTA code, written by Don Spong and +c modified by Jeremy Lore +c +c 7/2009 JL +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c PROGRAMMING NOTES +c -speed optimization +c -recursive ambipolar solver +c -fix numargs +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c 2024 Notes +c - Created PENTA_IMP_MAIN program and changed PENTA_IMP to a +c callable subroutine. +c +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + PROGRAM PENTA_MAIN + use penta_kind_mod + IMPLICIT NONE + integer(iknd) :: numargs, js, i_append + real(rknd) :: Er_min, Er_max, B_Eprl + character(100) :: arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8 + character(100) :: coeff_ext, run_ident + character(10) :: pprof_char + numargs = 8 + call getarg(1, arg1) + call getarg(2, arg2) + call getarg(3, arg3) + call getarg(4, arg4) + call getarg(5, arg5) + call getarg(6, arg6) + call getarg(7, arg7) + if (numargs .eq. 8) call getarg(8, arg8) + + !store command line args + coeff_ext = trim(adjustl(arg1)) + read(arg2,*) Er_min + read(arg3,*) Er_max + read(arg4,*) js + read(arg5,*) i_append + run_ident = trim(adjustl(arg6)) + pprof_char = trim(adjustl(arg7)) + if (numargs .eq. 8) read(arg8,*) B_Eprl + + ! Call subroutine + CALL PENTA_IMP(coeff_ext, Er_min, Er_max, js, i_append, + 1 run_ident, pprof_char, B_Eprl) + END PROGRAM PENTA_MAIN + From 8fbfe43a8290b1b0ca15c731205254a0c49caf73 Mon Sep 17 00:00:00 2001 From: lazersos Date: Thu, 22 Aug 2024 14:42:21 +0200 Subject: [PATCH 02/47] DKES: Added missing vimatrix dependency file. --- DKES/DKES.dep | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/DKES/DKES.dep b/DKES/DKES.dep index 11eb0160e..29eba5665 100644 --- a/DKES/DKES.dep +++ b/DKES/DKES.dep @@ -73,6 +73,9 @@ residue_dkes.o : \ vcmain.o : \ ../../LIBSTELL/$(LOCTYPE)/stel_kinds.o + + +vimatrix.o : vnamecl2.o : \ @@ -88,7 +91,7 @@ wrout.o : \ dkes_input.o \ dkes_realspace.o \ vnamecl2.o - + dkes_input_prepare.o : \ ../../LIBSTELL/$(LOCTYPE)/date_and_computer.o \ From 31a03cfb22ee9b1096d3179ecb385ff6fc12136c Mon Sep 17 00:00:00 2001 From: lazersos Date: Thu, 22 Aug 2024 14:43:20 +0200 Subject: [PATCH 03/47] THRIFT: Added interface to DKES (for testing). --- SHARE/make_all.inc | 10 + THRIFT/ObjectList | 1 + THRIFT/Sources/thrift_dkes.f90 | 308 ++++++++++++++++++++++++++++ THRIFT/Sources/thrift_input_mod.f90 | 11 +- THRIFT/Sources/thrift_vars.f90 | 8 + THRIFT/THRIFT.dep | 20 +- THRIFT/makethrift | 6 +- 7 files changed, 360 insertions(+), 4 deletions(-) create mode 100644 THRIFT/Sources/thrift_dkes.f90 diff --git a/SHARE/make_all.inc b/SHARE/make_all.inc index cb5dda4e8..730cf635c 100644 --- a/SHARE/make_all.inc +++ b/SHARE/make_all.inc @@ -111,6 +111,16 @@ else ifeq ($(MAKECMDGOALS),xthrift) $(BOOTSJ_DIR)/$(LOCTYPE)/$(LIB_BOOTSJ) \ $(BOOZ_DIR)/$(LOCTYPE)/$(LIB_BOOZ) \ $(DIAGNO_DIR)/$(LOCTYPE)/$(LIB_DIAGNO) + ifeq ($(LDKES),T) + DKES_DIR = ../../DKES + LIB_DKES = libdkes$(STATIC_EXT) + MOD_PATH += -I$(DKES_DIR)/$(LOCTYPE) + PRECOMP += -DDKES_OPT + LIB_LINK += $(DKES_DIR)/$(LOCTYPE)/$(LIB_DKES) + else + DKES_DIR = + LIB_DKES = + endif LGENE = F LCOILOPT = F LSFINCS = F diff --git a/THRIFT/ObjectList b/THRIFT/ObjectList index 983318595..634831d20 100644 --- a/THRIFT/ObjectList +++ b/THRIFT/ObjectList @@ -1,4 +1,5 @@ ObjectFiles = \ +thrift_dkes.o \ thrift_diagno.o \ thrift_run_diagnostics.o \ thrift_run_nbcd.o \ diff --git a/THRIFT/Sources/thrift_dkes.f90 b/THRIFT/Sources/thrift_dkes.f90 new file mode 100644 index 000000000..f55a45fb4 --- /dev/null +++ b/THRIFT/Sources/thrift_dkes.f90 @@ -0,0 +1,308 @@ +!----------------------------------------------------------------------- +! Subroutine: thrift_dkes +! Authors: S. Lazerson (samuel.lazerson@gauss-fusion.com) +! Date: 08/22/2024 +! Description: This subroutine calculates the DKES coefficients. +! see stellopt_dkes.f90 for more information. +!----------------------------------------------------------------------- + SUBROUTINE thrift_dkes(lscreen,iflag) +!----------------------------------------------------------------------- +! Libraries +!----------------------------------------------------------------------- + USE thrift_runtime + USE thrift_input_mod + USE thrift_vars, nrho_thrift => nrho + USE thrift_profiles_mod + USE thrift_equil + USE thrift_funcs + USE mpi_params + USE mpi_inc + USE read_wout_mod, ONLY: read_wout_file, read_wout_deallocate + USE read_boozer_mod, ONLY: bcast_boozer_vars + ! DKES LIBRARIES +#if defined(DKES_OPT) + USE Vimatrix + USE Vnamecl2 + USE Vcmain + USE dkes_input, lscreen_dkes => lscreen + USE dkes_realspace +#endif + +!----------------------------------------------------------------------- +! Subroutine Parameters +! iflag Error flag +!---------------------------------------------------------------------- + IMPLICIT NONE + LOGICAL, INTENT(in) :: lscreen + INTEGER, INTENT(inout) :: iflag +!----------------------------------------------------------------------- +! Local Variables +! ier Error flag +! iunit File unit number +!---------------------------------------------------------------------- + INTEGER :: i, j, k, l, istat, neqs, ier_phi, mystart, myend + INTEGER, DIMENSION(:), ALLOCATABLE :: ik_dkes + REAL(rprec), DIMENSION(:), ALLOCATABLE :: Earr_dkes, nuarr_dkes + REAL(rprec), DIMENSION(:), POINTER :: f0p1, f0p2, f0m1, f0m2 + REAL(rprec) :: tcpu0, tcpu1, tcpui, tcput, tcpu, tcpua, dkes_efield,& + phi_temp + CHARACTER :: tb*1 + CHARACTER*50 :: arg1(6) + INTEGER :: numargs, iodata, iout_opt, idata, iout + INTEGER :: m, n ! nmax, mmax (revmoed due to conflict with dkes_realspace + CHARACTER :: output_file*64, opt_file*64, dkes_input_file*64, temp_str*64 + LOGICAL :: lfirst_pass + +!----------------------------------------------- +! E x t e r n a l F u n c t i o n s +!----------------------------------------------- +!---------------------------------------------------------------------- +! BEGIN SUBROUTINE +!---------------------------------------------------------------------- + IF (iflag < 0) RETURN + IF (lscreen) WRITE(6,'(a)') ' -------------------- NEOCLASSICAL CALCULATION USING DKES -------------------' + IF (lvmec) THEN + ierr_mpi = 0 + ! First we deallocate the old arrays and allocate new ones +#if defined(MPI_OPT) + CALL MPI_BCAST(nruns_dkes,1,MPI_INTEGER,master,MPI_COMM_MYWORLD,ierr_mpi) +#endif + IF (ALLOCATED(DKES_L11p)) DEALLOCATE(DKES_L11p) + IF (ALLOCATED(DKES_L33p)) DEALLOCATE(DKES_L33p) + IF (ALLOCATED(DKES_L31p)) DEALLOCATE(DKES_L31p) + IF (ALLOCATED(DKES_L11m)) DEALLOCATE(DKES_L11m) + IF (ALLOCATED(DKES_L33m)) DEALLOCATE(DKES_L33m) + IF (ALLOCATED(DKES_L31m)) DEALLOCATE(DKES_L31m) + IF (ALLOCATED(DKES_scal11)) DEALLOCATE(DKES_scal11) + IF (ALLOCATED(DKES_scal33)) DEALLOCATE(DKES_scal33) + IF (ALLOCATED(DKES_scal31)) DEALLOCATE(DKES_scal31) + ALLOCATE(DKES_L11p(nruns_dkes),DKES_L33p(nruns_dkes),DKES_L31p(nruns_dkes),& + DKES_L11m(nruns_dkes),DKES_L33m(nruns_dkes),DKES_L31m(nruns_dkes),& + DKES_scal11(nruns_dkes),DKES_scal33(nruns_dkes),DKES_scal31(nruns_dkes)) + DKES_L11p=0.0; DKES_L33p=0.0; DKES_L31p=0.0 + DKES_L11m=0.0; DKES_L33m=0.0; DKES_L31m=0.0 + DKES_scal11=0.0; DKES_scal33=0.0; DKES_scal31=0.0 + l = 1 + ! Allocate Helper arrays to loop over dkes runs + ALLOCATE(ik_dkes(nruns_dkes),Earr_dkes(nruns_dkes),nuarr_dkes(nruns_dkes)) + IF (myworkid == master) THEN + DO k = 1, DKES_NS_MAX + DO i = 1, DKES_NSTAR_MAX + DO j = 1, DKES_NSTAR_MAX + IF ((DKES_K(k) > 0) .AND. (DKES_Erstar(i)<1E10) .AND. (DKES_Nustar(j) < 1E10)) THEN + ik_dkes(l) = DKES_K(k) + Earr_dkes(l) = DKES_Erstar(i) + nuarr_dkes(l) = DKES_Nustar(j) + l = l + 1 + END IF + END DO + END DO + END DO + ENDIF +#if defined(MPI_OPT) + ierr_mpi = 0 + CALL MPI_BCAST(ik_dkes,nruns_dkes,MPI_INTEGER,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(nuarr_dkes,nruns_dkes,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(Earr_dkes,nruns_dkes,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) +#endif + ! Read the wout file + CALL read_wout_file(proc_string, ier) + ! Now bcast the boozer parameters + CALL bcast_boozer_vars(master, MPI_COMM_MYWORLD, ierr_mpi) + ! Breakup work + CALL MPI_CALC_MYRANGE(MPI_COMM_MYWORLD,1,nruns_dkes,mystart,myend) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!! Parallel Work block + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Note there are multile file write statements which can + ! probably be removed so we only write what we need for PENTA + DO k = mystart,myend + ! Setup input + arg1(1) = TRIM(proc_string) + WRITE(arg1(2),'(i3)') ik_dkes(k) + WRITE(arg1(3),'(e20.10)') nuarr_dkes(k) + WRITE(arg1(4),'(e20.10)') dkes_efield + arg1(5) = 'F' + IF (lscreen .and. lfirst_pass) arg1(5) = 'T' + WRITE(temp_str,'(i3.3)') k + arg1(6) = '_s' // TRIM(temp_str) + ier_phi = 0 ! We don't read the boozmn or wout file we've done that already + CALL dkes_input_prepare(arg1,6,dkes_input_file,ier_phi) + output_file= 'dkesout.' // TRIM(proc_string) // '_s' // TRIM(temp_str) + opt_file= 'opt_dkes.' // TRIM(proc_string) // '_s' // TRIM(temp_str) !DAS 2/21/2000 !Probably won't need + summary_file = 'results.' // TRIM(proc_string) //'_s' // TRIM(temp_str) !record file addition + ! OPEN INPUT AND OUTPUT FILES FOR READING (AND WRITING OUTPUT) + idata = 7 + iout = 30 + iout_opt = 14 + iodata = idata + CALL safe_open(iodata, istat, dkes_input_file, 'old', 'formatted') + IF (istat .ne. 0) STOP 'Error reading input file in DKES' + ioout = iout + CALL safe_open(ioout, istat, output_file, 'replace', 'formatted') + IF (istat .ne. 0) STOP 'Error writing output file' + ioout_opt = iout_opt + CALL safe_open(ioout_opt, istat, opt_file, 'replace','formatted') + IF (istat .ne. 0) STOP 'Error writing opt_output file' + ! Read namelist (datain) input + lscreen_dkes = lscreen + nvalsb(1) = -bigint-1 + idisk = 1 + lfout = 0 + ibbi = 1 + borbi = 0 + READ (iodata, nml=dkes_indata, iostat=istat) + IF (istat .ne. 0) STOP 'Error reading dkes_indata NAMELIST in DKES' + CLOSE (iodata) + ! Recompute ntorb, mpolb for new style input where + ! borbi is input with actual index value, borbi(n,m) + IF (nvalsb(1) <= -bigint) THEN + nmax = -bigint; mmax = -bigint + DO n = -ntorbd, ntorbd + DO m = 0, mpolbd + IF (borbi(n,m) .ne. zero) THEN + nmax = max (nmax, abs(n)) + mmax = max (mmax, abs(m)) + END IF + END DO + END DO + ! User MAY input smaller values if he does not want to use entire array + ntorb = min (abs(ntorb), nmax) + mpolb = min (abs(mpolb-1), mmax) + ELSE + IF (ntorb > ntorbd) THEN + WRITE (ioout, 45) ntorb, ntorbd + STOP ' ntorb > ntorbd in DKES input' + END IF + IF (mpolb < 2) THEN + WRITE (ioout, 20) mpolb + STOP ' mpolb < 2 in DKES input' + ENDIF + END IF + IF (nzperiod .le. 0) nzperiod = 1 + IF (ipmb<0 .or. ipmb>2) ipmb = 0 + meshtz = MAX(0,meshtz) + !!!! END read_dkes_input + CALL set_mndim + neqs = mpnt*(lalpha + 1) + ALLOCATE (fzerop(2*neqs), fzerom(2*neqs), stat=istat) + !IF (istat .ne. 0) STOP 'allocation error(1) in dkes2!' + f0p1 => fzerop(1:neqs); f0p2 => fzerop(neqs+1:) !+ functions + f0m1 => fzerom(1:neqs); f0m2 => fzerom(neqs+1:) !- functions + CALL ftconv + CALL lcalc + CALL second0 (tcpu1); tcpui = tcpu1 - tcpu0 + CALL header + CALL second0 (tcpu0); tcput = zero + ! Here things get a bit screwy + ! DKES allows for an array of nu/v and E/v to be evaluated + ! but dkes_input_prepare does not. So nrun will always be 1 for stellopt + ! which is probably fine since we can parallelize over runs easier. + ! But in the future this logic could be improved. + nrun = MAX(nrun, 1) + nrun = MIN(nrun, krun) + DO i = 1, nrun + efield1 = efield(i) + cmul1 = cmul(i) + if(i .eq. 1) then + call safe_open(itab_out, istat, summary_file,'unknown', & + 'formatted') + write(itab_out,'("*",/,"cmul",a1,"efield",a1,"weov",a1,"wtov", & + & a1,"L11m",a1,"L11p",a1,"L31m",a1,"L31p",a1,"L33m",a1,"L33p", & + & a1,"scal11",a1,"scal13",a1,"scal33",a1,"max_residual", & + & a1,"chip",a1,"psip",a1,"btheta",a1,"bzeta",a1,"vp")') & + tb,tb,tb,tb,tb,tb,tb,tb,tb,tb,tb,tb,tb,tb,tb,tb,tb,tb + else if(i .gt. 1) then + open(itab_out,file=summary_file,status='unknown',position='append',form='formatted') + endif + CALL cescale (srces0) + WRITE (ioout, 950) dashes, cmul1, efield1, weov, wtov, wcyclo, vthermi + IF (ipmb < 2) THEN + iswpm = 1 + CALL blk5d (blk1, blk2, blk3, blk4, blk5, blk6, blk7, f0p1, f0p2, srces0) + IF (ier .ne. 0) THEN + WRITE (ioout, 1000) ier + STOP + ENDIF + CALL residue_dkes (blk1, blk2, blk3, blk4, blk5, blk6,& + f0p1, f0p2, srces0, rsd1p, rsd3p, g11p, g33p,& + g31p, g13p, crs1p, crs3p) + ENDIF + IF (ipmb .ne. 1) THEN + iswpm = 2 + CALL blk5d (blk1, blk2, blk3, blk4, blk5, blk6, blk7, f0m1, f0m2, srces0) + IF (ier .ne. 0) THEN + WRITE (ioout, 1050) ier + STOP + ENDIF + CALL residue_dkes (blk1, blk2, blk3, blk4, blk5, blk6,& + f0m1, f0m2, srces0, rsd1m, rsd3m, g11m, g33m,& + g31m, g13m, crs1m, crs3m) + ENDIF + ! This is a trick to get the arrays corretly sorted + DKES_rad_dex = i + IF (.not. lfirst_pass) lscreen_dkes = .FALSE. + CALL dkes_printout (f0p1, f0m1, f0p2, f0m2, srces0) + DKES_rad_dex = ik_dkes(i) + ! End trick + CALL second0 (tcpu1); tcpu = tcpu1 - tcpu0; tcpu0 = tcpu1; tcput = tcput + tcpu; tcpua = tcput/i + WRITE (ioout, 1100) tcpu + CLOSE(unit=itab_out) + END DO + !IF (lfout .ne. 0) CALL wrout (f0p1, f0m1, f0p2, f0m2, srces0) ! We don't need to do this + CALL free_mndim + DEALLOCATE (cols, al1, al2, al3, al4, bl1, bl2, bl3, bl4, cl1,& + cl2, cl3, cl4, cols0, omgl, al01, al02, al03, al04, bl01,& + bl02, bl03, bl04, cl01, cl02, cl03, cl04, fzerop, fzerom) + CLOSE(unit=ioout) + CLOSE(unit=ioout_opt) + lfirst_pass = .FALSE. + END DO + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!! Parallel Work block Done + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +#if defined(MPI_OPT) + CALL MPI_BARRIER(MPI_COMM_MYWORLD,ierr_mpi) + IF (myworkid == master) THEN + CALL MPI_REDUCE(MPI_IN_PLACE, DKES_L11p, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + CALL MPI_REDUCE(MPI_IN_PLACE, DKES_L11m, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + CALL MPI_REDUCE(MPI_IN_PLACE, DKES_L33p, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + CALL MPI_REDUCE(MPI_IN_PLACE, DKES_L33m, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + CALL MPI_REDUCE(MPI_IN_PLACE, DKES_L31p, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + CALL MPI_REDUCE(MPI_IN_PLACE, DKES_L31m, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + CALL MPI_REDUCE(MPI_IN_PLACE, DKES_scal11, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + CALL MPI_REDUCE(MPI_IN_PLACE, DKES_scal33, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + CALL MPI_REDUCE(MPI_IN_PLACE, DKES_scal31, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + ELSE + CALL MPI_REDUCE(DKES_L11p, DKES_L11p, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + CALL MPI_REDUCE(DKES_L11m, DKES_L11m, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + CALL MPI_REDUCE(DKES_L33p, DKES_L33p, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + CALL MPI_REDUCE(DKES_L33m, DKES_L33m, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + CALL MPI_REDUCE(DKES_L31p, DKES_L31p, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + CALL MPI_REDUCE(DKES_L31m, DKES_L31m, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + CALL MPI_REDUCE(DKES_scal11, DKES_scal11, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + CALL MPI_REDUCE(DKES_scal33, DKES_scal33, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + CALL MPI_REDUCE(DKES_scal31, DKES_scal31, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) + END IF +#endif + IF (lscreen) WRITE(6,'(a)') ' ----------------- DKES CALCULATION (DONE) ----------------' + IF (myworkid .ne. master) CALL read_wout_deallocate + + END IF + IF (lscreen) WRITE(6,'(a)') ' ------------------- DKES NEOCLASSICAL CALCULATION DONE ---------------------' + RETURN + 90 format(5e16.8) + 950 FORMAT(/9x,'CMUL',7x,'EFIELD',4x,'OMEGA-E/v',4x,'OMEGA-T/v',& + 5x,'OMEGA-ci',3x,'VI-THERMAL'/& + 7x,'(nu/v)',7x,'(Es/v)',2x,'(ExB drift)',4x,'(transit)',& + 3x, '(H+, B=B00)',2x,'(H+, Ti=1keV)'/,131a,/& + 1p,6e13.4/) + 20 FORMAT(' mpolb = ',i5,' is less than 2') + 45 FORMAT(' ntorb = ',i5,' is greater than ntorbd = ',i5) + 1000 FORMAT(/' blk5d error in block = ',i5) + 1050 FORMAT(/' blk5d error : ierm = ',i5) + 1100 FORMAT(/' time used: tcpu =',1p,e10.2,' sec'/) +!---------------------------------------------------------------------- +! END SUBROUTINE +!---------------------------------------------------------------------- + END SUBROUTINE thrift_dkes diff --git a/THRIFT/Sources/thrift_input_mod.f90 b/THRIFT/Sources/thrift_input_mod.f90 index 05e4862a5..d52d86b2d 100644 --- a/THRIFT/Sources/thrift_input_mod.f90 +++ b/THRIFT/Sources/thrift_input_mod.f90 @@ -35,7 +35,8 @@ MODULE thrift_input_mod targetposition_ecrh, rbeam_ecrh, & rfocus_ecrh, nra_ecrh, nphi_ecrh, & freq_ecrh, power_ecrh, & - pecrh_aux_t, pecrh_aux_f, ecrh_rc, ecrh_w + pecrh_aux_t, pecrh_aux_f, ecrh_rc, ecrh_w, & + dkes_k, dkes_Erstar, dkes_Nustar !----------------------------------------------------------------------- ! Subroutines @@ -84,6 +85,10 @@ SUBROUTINE init_thrift_input rfocus_ecrh = 0 nra_ecrh = 0 nphi_ecrh = 8 + ! DKES Vars + dkes_k = -1 + dkes_Erstar = 1E10 + dkes_Nustar = 1E10 RETURN END SUBROUTINE init_thrift_input @@ -115,7 +120,8 @@ SUBROUTINE read_thrift_input(filename, istat) CALL tolower(bootstrap_type) CALL tolower(eccd_type) leccd = eccd_type .ne. '' - nsj = nrho; + nsj = nrho + nruns_dkes = COUNT(dkes_k>0)*COUNT(dkes_Erstar<1E10)*COUNT(dkes_Nustar<1E10) RETURN END SUBROUTINE read_thrift_input @@ -150,6 +156,7 @@ SUBROUTINE write_thrift_namelist(iunit_out, istat) WRITE(iunit_out,outint) 'MBOZ',nboz WRITE(iunit_out,'(A)') '!---------- ECCD PARAMETERS ------------' WRITE(iunit_out,outstr) 'ECCD_TYPE',eccd_type + WRITE(iunit_out,'(A)') '!---------- DKES PARAMETERS ------------' WRITE(iunit_out,'(A)') '/' RETURN END SUBROUTINE write_thrift_namelist diff --git a/THRIFT/Sources/thrift_vars.f90 b/THRIFT/Sources/thrift_vars.f90 index c6d820fc1..21290a09c 100644 --- a/THRIFT/Sources/thrift_vars.f90 +++ b/THRIFT/Sources/thrift_vars.f90 @@ -110,5 +110,13 @@ MODULE thrift_vars REAL(rprec), DIMENSION(nsys,3) :: antennaposition_ecrh, & targetposition_ecrh,rbeam_ecrh,rfocus_ecrh + ! For DKES + INTEGER, PARAMETER :: DKES_NS_MAX = 64 + INTEGER, PARAMETER :: DKES_NSTAR_MAX = 32 + INTEGER :: nruns_dkes + INTEGER, DIMENSION(:), POINTER :: DKES_rundex + INTEGER, DIMENSION(DKES_NS_MAX) :: DKES_K + REAL(rprec), DIMENSION(DKES_NSTAR_MAX) :: DKES_Erstar, DKES_Nustar + END MODULE thrift_vars diff --git a/THRIFT/THRIFT.dep b/THRIFT/THRIFT.dep index f93470d1c..38ac7e9a9 100644 --- a/THRIFT/THRIFT.dep +++ b/THRIFT/THRIFT.dep @@ -22,7 +22,25 @@ thrift_boozer.o : \ $(BOOZ_DIR)/$(LOCTYPE)/booz_persistent.o \ $(LIB_DIR)/$(LOCTYPE)/read_boozer_mod.o \ $(LIB_DIR)/$(LOCTYPE)/mpi_inc.o \ - $(LIB_DIR)/$(LOCTYPE)/mpi_params.o + $(LIB_DIR)/$(LOCTYPE)/mpi_params.o + +thrift_dkes.o : \ + thrift_runtime.o \ + thrift_vars.o \ + thrift_input_mod.o \ + thrift_profiles_mod.o \ + thrift_equil.o \ + thrift_funcs.o \ + $(LIB_DIR)/$(LOCTYPE)/mpi_inc.o \ + $(LIB_DIR)/$(LOCTYPE)/mpi_params.o \ + $(LIB_DIR)/$(LOCTYPE)/bootsj_input.o \ + $(DKES_DIR)/$(LOCTYPE)/vimatrix.o \ + $(DKES_DIR)/$(LOCTYPE)/vnamecl2.o \ + $(DKES_DIR)/$(LOCTYPE)/vcmain.o \ + $(DKES_DIR)/$(LOCTYPE)/dkes_input.o \ + $(DKES_DIR)/$(LOCTYPE)/dkes_realspace.o \ + $(LIB_DIR)/$(LOCTYPE)/read_wout_mod.o \ + $(LIB_DIR)/$(LOCTYPE)/read_boozer_mod.o thrift_diagno.o : \ thrift_runtime.o \ diff --git a/THRIFT/makethrift b/THRIFT/makethrift index 0e3442056..fd2ad7315 100644 --- a/THRIFT/makethrift +++ b/THRIFT/makethrift @@ -15,7 +15,7 @@ VPATH = $(SPATH) .SUFFIXES : .SUFFIXES : .f .f90 .o -xthrift: $(LIB) $(LIB_VMEC) $(LIB_BOOTSJ) $(LIB_BOOZ) $(LIB_DIAGNO) $(ObjectFiles) +xthrift: $(LIB) $(LIB_VMEC) $(LIB_BOOTSJ) $(LIB_BOOZ) $(LIB_DIAGNO) $(LIB_DKES) $(ObjectFiles) $(LINK) $@ $(ObjectFiles) $(LIB_LINK) ifdef VMEC_DIR # @rm $(VMEC_DIR)/$(LOCTYPE)/$(LIB_VMEC) @@ -89,3 +89,7 @@ $(LIB_BOOZ) : #Construct diagno library. $(LIB_DIAGNO) : @cd $(DIAGNO_DIR); make $(TYPE); cd $(LOCTYPE); ar -cruv $(LIB_DIAGNO) *.o + +#Construct DKES library. +$(LIB_DKES) : + @cd $(DKES_DIR); make $(TYPE); cd $(LOCTYPE); ar -cruv $(LIB_DKES) *.o From fcd45a15e92025fe0e4a90b27630edffb464e6a5 Mon Sep 17 00:00:00 2001 From: lazersos Date: Sat, 24 Aug 2024 10:45:15 +0900 Subject: [PATCH 04/47] PENTA: Make possible to read namelist from input.ext file. --- PENTA/Sources/penta_imp.f | 50 ++++++++++++++++++++++++++++++++--- PENTA/Sources/penta_modules.f | 3 ++- 2 files changed, 49 insertions(+), 4 deletions(-) diff --git a/PENTA/Sources/penta_imp.f b/PENTA/Sources/penta_imp.f index b4acb8633..b37fb7632 100644 --- a/PENTA/Sources/penta_imp.f +++ b/PENTA/Sources/penta_imp.f @@ -9,6 +9,7 @@ SUBROUTINE penta_imp(coeff_ext, Er_min, Er_max, js, i_append, use Er_roots_pass use io_unit_spec use parameter_pass + USE safe_open_mod implicit none logical, parameter :: input_is_Er = .true. !If true, Er range is (V/cm) otherwise eEr/kT_e @@ -33,7 +34,7 @@ SUBROUTINE penta_imp(coeff_ext, Er_min, Er_max, js, i_append, integer(iknd) :: num_species, iroot, 1 num_ion_species, inv_err, nn_inv, j, ispec, tmp_ind, ie, - 2 ispec1, ispec2, spec1_ind, spec2_ind, numargs + 2 ispec1, ispec2, spec1_ind, spec2_ind, numargs, istat real(rknd) :: vth_e, loglambda, sigma_S, Er_test, U2, 1 lX_sum1, lX_sum2, J_E_tot, J_E_cl, B_Eprl, @@ -57,6 +58,8 @@ SUBROUTINE penta_imp(coeff_ext, Er_min, Er_max, js, i_append, character(100) :: arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, 1 arg9, format_tmp, str_num, fstatus, fpos + LOGICAL :: lexist + namelist / ion_params / num_ion_species,Z_ion_init,miomp_init !external functions @@ -80,9 +83,50 @@ SUBROUTINE penta_imp(coeff_ext, Er_min, Er_max, js, i_append, IF (PRESENT(B_Eprl_in)) B_Eprl= B_Eprl_in !Read namelist file for ion parameters - open(iu_nl,file="ion_params",status="old",form="formatted") - read(iu_nl,nml=ion_params) + INQUIRE(FILE='input.'//TRIM(run_ident),EXIST=lexist) + IF (lexist) THEN + CALL safe_open(iu_nl,istat,'input.'//TRIM(run_ident),'old', + 1 'formatted') + ELSE + CALL safe_open(iu_nl,istat,'ion_params','old','formatted') + WRITE(6,*) + 1 ' -Reading ion_params namelist from ion_params file' + END IF + READ(iu_nl,NML=ion_params,IOSTAT=istat) + IF (istat /= 0) THEN + CLOSE(iu_nl) +! IF (lverb) THEN +! WRITE(6,*) +! 1 ' -Could not find ion_params namelist in input.'// +! 2 TRIM(run_ident)//' file' +! WRITE(6,*) +! 1 ' -Attempting to read ion_params namelist from ion_params file' +! END IF + iu_nl = 25 + istat = 0 + CALL safe_open(iu_nl,istat,'ion_params','old','formatted') + IF (istat /= 0) THEN + WRITE(6,*) ' -Could not open ion_params file' + !WRITE(6,*) ' -Dumping ion_params to screen' + stop + END IF + istat = 0 + READ(iu_nl,NML=ion_params,IOSTAT=istat) + IF (istat /= 0) THEN + WRITE(6,*) + 1 ' -Could not read ion_params from ion_params file' + !WRITE(6,*) ' -Dumping ion_params to screen' + !CALL write_diagno_input(6,istat) + !backspace(iu_nl) + !read(iu_nl,fmt='(A)') line + !write(6,'(A)') 'Invalid line in namelist: '//TRIM(line) + stop + END IF close(iu_nl) + END IF + !open(iu_nl,file="ion_params",status="old",form="formatted") + !read(iu_nl,nml=ion_params) + !close(iu_nl) !Get command line arguments !numargs = iargc() diff --git a/PENTA/Sources/penta_modules.f b/PENTA/Sources/penta_modules.f index b0daba464..f38b958f7 100644 --- a/PENTA/Sources/penta_modules.f +++ b/PENTA/Sources/penta_modules.f @@ -27,7 +27,8 @@ end module phys_const c----------------------------------------------------------------- c module io_unit_spec - integer, parameter :: iu_nl=21, iu_vmec=22, iu_pprof=23, + INTEGER :: iu_nl = 21 + integer, parameter :: iu_vmec=22, iu_pprof=23, 1 iu_coeff=24, iu_Ufile=25, iu_flux_out=10, iu_pprof_out=11, 2 iu_fvEr_out=12, iu_flows_out=13 end module io_unit_spec From 3f2aff746bef2fa675cb3abe28efefe4949734d7 Mon Sep 17 00:00:00 2001 From: lazersos Date: Sat, 24 Aug 2024 10:45:45 +0900 Subject: [PATCH 05/47] THRIFT: Now links to PENTA sources --- SHARE/make_all.inc | 8 ++++++-- THRIFT/makethrift | 9 ++++++++- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/SHARE/make_all.inc b/SHARE/make_all.inc index 730cf635c..cd26367da 100644 --- a/SHARE/make_all.inc +++ b/SHARE/make_all.inc @@ -72,7 +72,7 @@ ifeq ($(MAKECMDGOALS),xstelloptv2) $(DIAGNO_DIR)/$(LOCTYPE)/$(LIB_DIAGNO) \ $(FIELDLINES_DIR)/$(LOCTYPE)/$(LIB_FIELDLINES) \ $(JINV_DIR)/$(LOCTYPE)/$(LIB_JINV) \ - $(MGRID_DIR)/$(LOCTYPE)/$(LIB_MGRID) + $(MGRID_DIR)/$(LOCTYPE)/$(LIB_MGRID) ifeq ($(LDKES),T) DKES_DIR = ../../DKES LIB_DKES = libdkes$(STATIC_EXT) @@ -107,10 +107,14 @@ else ifeq ($(MAKECMDGOALS),xthrift) DIAGNO_DIR = ../../DIAGNO LIB_DIAGNO= libdiagno$(STATIC_EXT) MOD_PATH+= -I$(DIAGNO_DIR)/$(LOCTYPE) + PENTA_DIR= ../../PENTA + LIB_PENTA= libpenta$(STATIC_EXT) + MOD_PATH+= -I$(PENTA_DIR)/$(LOCTYPE) LIB_LINK= $(VMEC_DIR)/$(LOCTYPE)/$(LIB_VMEC) \ $(BOOTSJ_DIR)/$(LOCTYPE)/$(LIB_BOOTSJ) \ $(BOOZ_DIR)/$(LOCTYPE)/$(LIB_BOOZ) \ - $(DIAGNO_DIR)/$(LOCTYPE)/$(LIB_DIAGNO) + $(DIAGNO_DIR)/$(LOCTYPE)/$(LIB_DIAGNO) \ + $(PENTA_DIR)/$(LOCTYPE)/$(LIB_PENTA) ifeq ($(LDKES),T) DKES_DIR = ../../DKES LIB_DKES = libdkes$(STATIC_EXT) diff --git a/THRIFT/makethrift b/THRIFT/makethrift index fd2ad7315..b221d5249 100644 --- a/THRIFT/makethrift +++ b/THRIFT/makethrift @@ -15,7 +15,7 @@ VPATH = $(SPATH) .SUFFIXES : .SUFFIXES : .f .f90 .o -xthrift: $(LIB) $(LIB_VMEC) $(LIB_BOOTSJ) $(LIB_BOOZ) $(LIB_DIAGNO) $(LIB_DKES) $(ObjectFiles) +xthrift: $(LIB) $(LIB_VMEC) $(LIB_BOOTSJ) $(LIB_BOOZ) $(LIB_DIAGNO) $(LIB_DKES) $(LIB_PENTA) $(ObjectFiles) $(LINK) $@ $(ObjectFiles) $(LIB_LINK) ifdef VMEC_DIR # @rm $(VMEC_DIR)/$(LOCTYPE)/$(LIB_VMEC) @@ -29,6 +29,9 @@ endif ifdef DIAGNO_DIR @rm $(DIAGNO_DIR)/$(LOCTYPE)/$(LIB_DIAGNO) endif +ifdef PENTA_DIR + @rm $(PENTA_DIR)/$(LOCTYPE)/$(LIB_PENTA) +endif #Compile object files defined in OBJLIST. .f.o : @@ -93,3 +96,7 @@ $(LIB_DIAGNO) : #Construct DKES library. $(LIB_DKES) : @cd $(DKES_DIR); make $(TYPE); cd $(LOCTYPE); ar -cruv $(LIB_DKES) *.o + +#Construct PENTA library. +$(LIB_PENTA) : + @cd $(PENTA_DIR); make $(TYPE); cd $(LOCTYPE); ar -cruv $(LIB_PENTA) *.o From 96ef2082aaf0baa6e0d95c7379f4441cda0cbaac Mon Sep 17 00:00:00 2001 From: lazersos Date: Mon, 26 Aug 2024 10:08:35 +0900 Subject: [PATCH 06/47] PENTA: Fixed call to vmec_2 routine. --- PENTA/Sources/penta_imp.f | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/PENTA/Sources/penta_imp.f b/PENTA/Sources/penta_imp.f index b37fb7632..7553a7c5b 100644 --- a/PENTA/Sources/penta_imp.f +++ b/PENTA/Sources/penta_imp.f @@ -173,7 +173,8 @@ SUBROUTINE penta_imp(coeff_ext, Er_min, Er_max, js, i_append, allocate(u_theta(num_species),u_zeta(num_species)) !read input files - call read_vmec_file(js,run_ident) + !call read_vmec_file(js,run_ident) + call read_vmec_file_2(js,run_ident) call read_pprof_file(js,pprof_char,num_ion_species,roa_surf,arad) call read_lmn_star_files(coeff_ext) call read_Utilde2_file(roa_surf,U2) From e7aca022fa09f1f10273a8f588e196ce4181add0 Mon Sep 17 00:00:00 2001 From: lazersos Date: Mon, 26 Aug 2024 10:08:56 +0900 Subject: [PATCH 07/47] PENTA: Added routine to set profiles via subroutine call. --- PENTA/Sources/read_input_file_mod.f | 97 ++++++++++++++++++++++++++++- 1 file changed, 94 insertions(+), 3 deletions(-) diff --git a/PENTA/Sources/read_input_file_mod.f b/PENTA/Sources/read_input_file_mod.f index 3b00accf8..bf9509ad6 100644 --- a/PENTA/Sources/read_input_file_mod.f +++ b/PENTA/Sources/read_input_file_mod.f @@ -70,8 +70,8 @@ subroutine read_vmec_file(js,run_ident) integer(iknd) :: js character(60) :: run_ident !local variables - integer(iknd) :: j, js_min, js_max - integer(iknd), dimension(:), allocatable :: js_vmec + integer(iknd) :: j, js_min, js_max + integer(iknd), dimension(:), allocatable :: js_vmec real(rknd), dimension(:), allocatable :: r_vmec, roa_vmec 1 ,chip_vmec, psip_vmec, btheta_vmec, bzeta_vmec, vp_vmec, bsq_vmec 2 ,iota_vmec @@ -125,7 +125,7 @@ subroutine read_vmec_file_2(js,run_ident) character(60) :: run_ident !local variables integer(iknd) :: j, js_min, js_max, ierr, u, v, mn - integer(iknd), dimension(:), allocatable :: js_vmec + integer(iknd), dimension(:), allocatable :: js_vmec ! real(rknd), dimension(:), allocatable :: r_vmec, roa_vmec ! 1 ,chip_vmec, psip_vmec, btheta_vmec, bzeta_vmec, vp_vmec, bsq_vmec ! 2 ,iota_vmec @@ -284,6 +284,97 @@ subroutine read_pprof_file(js,pprof_char,nis,roa_surf,arad) end subroutine read_pprof_file c c---------------------------------------------------------------------------- +c Reads the file "plasma_profiles*.dat" and assigns plasma parameters for the +c current surface. Variables are passed using module pprof_pass +c---------------------------------------------------------------------------- +c + subroutine set_pprof(np_prof,nis,roa_prof,ne_prof,Te_prof, + 1 ni_local,Ti_local,roa_surf,arad) + use penta_kind_mod + use pprof_pass + use bspline + use io_unit_spec + implicit none + ! Input Parameters + INTEGER(iknd) :: np_prof + INTEGER(iknd) :: nis + REAL(rknd), DIMENSION(np_prof) :: roa_prof + REAL(rknd), DIMENSION(np_prof) :: ne_prof + REAL(rknd), DIMENSION(np_prof) :: Te_prof + REAL(rknd), DIMENSION(np_prof,nis) :: ni_local + REAL(rknd), DIMENSION(np_prof,nis) :: Ti_local + INTEGER(iknd) :: js + REAL(rknd) :: roa_surf + REAL(rknd) :: arad + integer(iknd) :: ispec, kord_prof + real(rknd), dimension(:), allocatable :: + 1 Te_knot_array, spl_Te, ne_knot_array, + 2 spl_ne, ni_knot_array, spl_ni, Ti_knot_array, spl_Ti + + ! Some global helpers for ions (Not sure these are really needed) + IF (ALLOCATED(ni_prof)) DEALLOCATE(ni_prof) + IF (ALLOCATED(Ti_prof)) DEALLOCATE(Ti_prof) + ALLOCATE(ni_prof(np_prof,nis),Ti_prof(np_prof,nis)) + ni_prof = ni_local + Ti_prof = Ti_local + + !spline fit ne,Te,Ti,ni profiles and evaluate at r/a of current surface + kord_prof = 3 !spline order for profile fitting + allocate(ne_knot_array(np_prof+kord_prof),spl_ne(np_prof)) + allocate(Te_knot_array(np_prof+kord_prof),spl_Te(np_prof)) + allocate(ni_knot_array(np_prof+kord_prof),spl_ni(np_prof)) + allocate(Ti_knot_array(np_prof+kord_prof),spl_Ti(np_prof)) + + !fit electron profiles + call dbsnak(np_prof,roa_prof,kord_prof,ne_knot_array) + call dbsnak(np_prof,roa_prof,kord_prof,Te_knot_array) + call dbsint(np_prof,roa_prof,ne_prof,kord_prof, + 1 ne_knot_array,spl_ne) + call dbsint(np_prof,roa_prof,te_prof,kord_prof, + 1 Te_knot_array,spl_Te) + + !loop over ion species, fit profiles and assign values + do ispec=1,nis + call dbsnak(np_prof,roa_prof,kord_prof,ti_knot_array) + call dbsnak(np_prof,roa_prof,kord_prof,ni_knot_array) + call dbsint(np_prof,roa_prof,ni_prof(:,ispec),kord_prof, + 1 ni_knot_array,spl_ni) + call dbsint(np_prof,roa_prof,Ti_prof(:,ispec),kord_prof, + 1 Ti_knot_array,spl_ti) + !evaluate spline fit at r/a of the test surface + ni(ispec) = dbsval(roa_surf,kord_prof,ni_knot_array,np_prof, + 1 spl_ni) + Ti(ispec) = dbsval(roa_surf,kord_prof,ti_knot_array,np_prof, + 1 spl_ti) + !evaluate derivatives (d/dr) at r/a + dnidr(ispec) = dbsder(1,roa_surf,kord_prof,ni_knot_array, + 1 np_prof,spl_ni)/arad + dTidr(ispec) = dbsder(1,roa_surf,kord_prof,Ti_knot_array, + 1 np_prof,spl_Ti)/arad + enddo + + !evaluate spline fit at r/a of the test surface for electrons + ne = dbsval(roa_surf,kord_prof,ne_knot_array,np_prof,spl_ne) + Te = dbsval(roa_surf,kord_prof,te_knot_array,np_prof,spl_te) + !evaluate derivatives (d/dr) at r/a for electrons + dnedr = dbsder(1,roa_surf,kord_prof,ne_knot_array, + 1 np_prof,spl_ne)/arad + dTedr = dbsder(1,roa_surf,kord_prof,Te_knot_array, + 1 np_prof,spl_Te)/arad + !convert units to mks + ne=ne*1.e18_rknd + ni=ni*1.e18_rknd + dnedr=dnedr*1.e18_rknd + dnidr=dnidr*1.e18_rknd + + !deallocate variables + deallocate(Te_knot_array,spl_Te) + deallocate(ne_knot_array,spl_ne) + deallocate(ni_knot_array,spl_ni) + deallocate(Ti_knot_array, spl_Ti) + end subroutine set_pprof +c +c---------------------------------------------------------------------------- c Reads the file *star_lijs_*** files and assigns variables for passing c using module coeff_var_pass c---------------------------------------------------------------------------- From 7cd76c88b47efea6b9e599e2fc4f235ae8ed5e27 Mon Sep 17 00:00:00 2001 From: lazersos Date: Wed, 28 Aug 2024 10:34:34 +0900 Subject: [PATCH 08/47] PENTA: Added routine to set LMNIJ coefficients. --- PENTA/Sources/read_input_file_mod.f | 48 +++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/PENTA/Sources/read_input_file_mod.f b/PENTA/Sources/read_input_file_mod.f index bf9509ad6..22e652943 100644 --- a/PENTA/Sources/read_input_file_mod.f +++ b/PENTA/Sources/read_input_file_mod.f @@ -465,5 +465,53 @@ subroutine read_lmn_star_files(coeff_ext) num_cm=size(cmul_ms); num_em=size(efield_ms) num_cn=size(cmul_ns); num_en=size(efield_ns) end subroutine read_lmn_star_files +c +c---------------------------------------------------------------------------- +c Assigns variables for passing using module coeff_var_pass +c---------------------------------------------------------------------------- +c + subroutine set_lmn_star_files(nc,ne,ifile,cmul_vec,efield_vec, + 1 coef2d) + use penta_kind_mod + use coeff_var_pass + use io_unit_spec + implicit none + ! Input parameters + INTEGER, INTENT(IN) :: nc + INTEGER, INTENT(IN) :: ne + INTEGER, INTENT(IN) :: ifile ! 1=L; 2=M; 3=N + REAL(rknd), DIMENSION(nc), INTENT(IN) :: cmul_vec + REAL(rknd), DIMENSION(ne), INTENT(IN) :: efield_vec + REAL(rknd), DIMENSION(nc,ne), INTENT(IN) :: coef2d + + IF (ifile == 1) THEN + allocate(cmul_ls(nc)) + allocate(efield_ls(ne)) + allocate(coef2d_ls(nc,ne)) + cmul_ls=cmul_vec + efield_ls=efield_vec + coef2d_ls=coef2d + num_cl=size(cmul_ls); num_el=size(efield_ls) + ELSEIF (ifile == 2) THEN + allocate(cmul_ms(nc)) + allocate(efield_ms(ne)) + allocate(coef2d_ms(nc,ne)) + cmul_ms=cmul_vec + efield_ms=efield_vec + coef2d_ms=coef2d + num_cm=size(cmul_ms); num_em=size(efield_ms) + ELSEIF (ifile == 3) THEN + allocate(cmul_ns(nc)) + allocate(efield_ns(ne)) + allocate(coef2d_ns(nc,ne)) + cmul_ns=cmul_vec + efield_ns=efield_vec + coef2d_ns=coef2d + num_cn=size(cmul_ns); num_en=size(efield_ns) + ELSE + WRITE(6,*) 'set_lmn_star_files: Unknown ifile=',ifile + END IF + RETURN + end subroutine set_lmn_star_files end module read_input_file_mod \ No newline at end of file From 9a3974457d5e5778b1448d7173a1a1bb45f7b13f Mon Sep 17 00:00:00 2001 From: lazersos Date: Wed, 28 Aug 2024 10:35:26 +0900 Subject: [PATCH 09/47] PENTA: Cleanup of routine a bit. --- PENTA/Sources/penta_imp.f | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/PENTA/Sources/penta_imp.f b/PENTA/Sources/penta_imp.f index 7553a7c5b..31477113b 100644 --- a/PENTA/Sources/penta_imp.f +++ b/PENTA/Sources/penta_imp.f @@ -339,7 +339,6 @@ SUBROUTINE penta_imp(coeff_ext, Er_min, Er_max, js, i_append, !write fluxes vs Er write(str_num,*) num_ion_species*2+3 - !write(str_num,'(i)') num_ion_species*2+3 format_tmp='(f7.4,' // trim(adjustl(str_num)) // '(" ",e15.7))' write(iu_fvEr_out,format_tmp) roa_surf,Er_test/100._rknd, 1 gamma_e(ie),gamma_i(:,ie),q_e(ie),q_i(:,ie) @@ -430,7 +429,6 @@ SUBROUTINE penta_imp(coeff_ext, Er_min, Er_max, js, i_append, !write fluxes to file "flux_vs_roa" write(str_num,*) 2*num_ion_species+5 - !write(str_num,'(i)') 2*num_ion_species+5 format_tmp='(f7.4,' // trim(adjustl(str_num)) // '(" ",e15.7))' write(iu_flux_out,format_tmp) roa_surf,Er_test/100._rknd, 1 eaEr_o_kTe,gamma_e_ambi(iroot),q_e_ambi(iroot), @@ -438,7 +436,6 @@ SUBROUTINE penta_imp(coeff_ext, Er_min, Er_max, js, i_append, !write flows to file "flows_vs_roa" write(str_num,*) 4*num_ion_species+6 - !write(str_num,'(i)') 4*num_ion_species+6 format_tmp='(f7.4,' // trim(adjustl(str_num)) // '(" ",e15.7))' write(iu_flows_out,format_tmp) roa_surf,Er_test/100._rknd, 1 eaEr_o_kTe,B_uprle(iroot),B_uprli(:,iroot),B_qprle(iroot), @@ -447,16 +444,13 @@ SUBROUTINE penta_imp(coeff_ext, Er_min, Er_max, js, i_append, !write output files - !"plasma_profiles_check" - !write(str_num,'(i)') 4*num_ion_species+5 - write(str_num,*) 4*num_ion_species+5 + write(str_num,*) 4*num_ion_species+5 format_tmp='(' // trim(adjustl(str_num)) // '(" ",e15.7))' write(iu_pprof_out,format_tmp) roa_surf,Te,ne,dnedr,dTedr, 1 Ti,ni,dnidr,dTidr !write to screen - !write(str_num,'(i)') num_roots - write(str_num,*) num_roots + write(str_num,*) num_roots format_tmp='(f7.3,' // trim(adjustl(str_num)) // '(" ",e15.4))' write(*,format_tmp) roa_surf,er_roots/100._rknd From 0c45d24cdf16182f5d5d47a6ac847af76a201fca Mon Sep 17 00:00:00 2001 From: lazersos Date: Wed, 28 Aug 2024 11:01:51 +0900 Subject: [PATCH 10/47] PENTA: Changed name for VMEC reading routine. --- PENTA/Sources/penta_imp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PENTA/Sources/penta_imp.f b/PENTA/Sources/penta_imp.f index 31477113b..b88025e37 100644 --- a/PENTA/Sources/penta_imp.f +++ b/PENTA/Sources/penta_imp.f @@ -174,7 +174,7 @@ SUBROUTINE penta_imp(coeff_ext, Er_min, Er_max, js, i_append, !read input files !call read_vmec_file(js,run_ident) - call read_vmec_file_2(js,run_ident) + call read_vmec_file_wout(js,run_ident) call read_pprof_file(js,pprof_char,num_ion_species,roa_surf,arad) call read_lmn_star_files(coeff_ext) call read_Utilde2_file(roa_surf,U2) From a4cb861007488e7ad016c23c5f9eaf750c8d82bd Mon Sep 17 00:00:00 2001 From: lazersos Date: Wed, 28 Aug 2024 11:02:29 +0900 Subject: [PATCH 11/47] PENTA: Changed VMEC reading routine name and added Utilde2 setting routine. --- PENTA/Sources/penta_main.f | 9 ++++ PENTA/Sources/read_input_file_mod.f | 83 ++++++++++++++++++----------- 2 files changed, 60 insertions(+), 32 deletions(-) diff --git a/PENTA/Sources/penta_main.f b/PENTA/Sources/penta_main.f index 541f6a2dd..e49aa383a 100644 --- a/PENTA/Sources/penta_main.f +++ b/PENTA/Sources/penta_main.f @@ -88,11 +88,16 @@ c 2024 Notes c - Created PENTA_IMP_MAIN program and changed PENTA_IMP to a c callable subroutine. +c - Can read wout files instead of profile_data_**** +c - Routine to set profiles instead of reading plasma_profiles_**** +c - Routine to set L/M/NIJ coefficients instead of reading them. +c - Routine to set Utilde2 value instead of reading from file. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc PROGRAM PENTA_MAIN use penta_kind_mod + USE read_wout_mod, ONLY: read_wout_deallocate IMPLICIT NONE integer(iknd) :: numargs, js, i_append real(rknd) :: Er_min, Er_max, B_Eprl @@ -122,5 +127,9 @@ PROGRAM PENTA_MAIN ! Call subroutine CALL PENTA_IMP(coeff_ext, Er_min, Er_max, js, i_append, 1 run_ident, pprof_char, B_Eprl) + + ! DEALLOCATE VMEC DATA + CALL read_wout_deallocate + END PROGRAM PENTA_MAIN diff --git a/PENTA/Sources/read_input_file_mod.f b/PENTA/Sources/read_input_file_mod.f index 22e652943..bc204c093 100644 --- a/PENTA/Sources/read_input_file_mod.f +++ b/PENTA/Sources/read_input_file_mod.f @@ -57,6 +57,40 @@ subroutine read_Utilde2_file(roa_test,U2) end subroutine read_Utilde2_file c c---------------------------------------------------------------------------- +c Sets the file "Utilde2_profile" and returns for the run surface +c---------------------------------------------------------------------------- +c + subroutine set_Utilde2_file(num_pts,roa_Utilde,U_tilde2, + 1 roa_test,U2) + use penta_kind_mod + use io_unit_spec + implicit none + ! Input parameters + INTEGER(iknd), INTENT(IN) :: num_pts + REAL(rknd), DIMENSION(num_pts), INTENT(IN) :: roa_Utilde + REAL(rknd), DIMENSION(num_pts), INTENT(IN) :: U_tilde2 + REAL(rknd), INTENT(IN) :: roa_test + REAL(rknd), INTENT(OUT) :: U2 + !dummy variables + !local varibles + integer(iknd) :: j, js_surf, num_pts, kord_prof + real(rknd), dimension(:), allocatable :: U2_knot_array, spl_U2 + + !spline fit profile and evaluate at r/a of current surface + kord_prof = 3 !spline order + allocate(U2_knot_array(num_pts+kord_prof)) + allocate(spl_U2(num_pts)) + call dbsnak(num_pts,roa_Utilde,kord_prof,U2_knot_array) + call dbsint(num_pts,roa_Utilde,U_tilde2,kord_prof, + 1 U2_knot_array,spl_U2) + U2 = dbsval(roa_test,kord_prof,U2_knot_array,num_pts, + 1 spl_U2) + + !deallocate variables + deallocate(U2_knot_array,spl_U2) + end subroutine set_Utilde2_file +c +c---------------------------------------------------------------------------- c Reads the file "profile_data_***" and assigns VMEC variables, passed c using module vmec_var_pass c---------------------------------------------------------------------------- @@ -108,7 +142,7 @@ end subroutine read_vmec_file c Reads the wout file for vmec information (SAL 05/30/13) c---------------------------------------------------------------------------- c - subroutine read_vmec_file_2(js,run_ident) + subroutine read_vmec_file_wout(js,run_ident) use vmec_var_pass use penta_kind_mod use io_unit_spec @@ -118,7 +152,8 @@ subroutine read_vmec_file_2(js,run_ident) 3 Rmajor_vmec => Rmajor, 4 Aminor_vmec => Aminor, 5 phi_vmec => phi, bmnc, mnmax_nyq, gmnc, - 6 xm_nyq, xn_nyq, ns_vmec => ns + 6 xm_nyq, xn_nyq, ns_vmec => ns, + 7 lwout_opened implicit none !dummy variables integer(iknd) :: js @@ -126,37 +161,23 @@ subroutine read_vmec_file_2(js,run_ident) !local variables integer(iknd) :: j, js_min, js_max, ierr, u, v, mn integer(iknd), dimension(:), allocatable :: js_vmec -! real(rknd), dimension(:), allocatable :: r_vmec, roa_vmec -! 1 ,chip_vmec, psip_vmec, btheta_vmec, bzeta_vmec, vp_vmec, bsq_vmec -! 2 ,iota_vmec REAL(rknd) :: TWOPI, top, bottom, theta, zeta, arg, cs_arg, 1 jacob_vmec, bbf character(60) :: vmec_fname, ch_dum, tb = char(9) + INTEGER(iknd), PARAMETER :: nu_int = 360 + INTEGER(iknd), PARAMETER :: nv_int = 360 + ! Initializations TWOPI = 8*ATAN(1.0_rknd) - !Read VMEC profile data file -! vmec_fname = "profile_data_" // run_ident -! open(unit=iu_vmec,file=vmec_fname,status="old") - CALL read_wout_file(TRIM(run_ident),ierr) - -! read(iu_vmec,*) js_min, js_max -! allocate(js_vmec(js_max),r_vmec(js_max),roa_vmec(js_max)) -! allocate(chip_vmec(js_max),psip_vmec(js_max),btheta_vmec(js_max)) -! allocate(bzeta_vmec(js_max),vp_vmec(js_max),bsq_vmec(js_max)) -! allocate(iota_vmec(js_max)) -! read(iu_vmec,*) arad, Rmajor -! read(iu_vmec,'(a10)') ch_dum -! do j = js_min,js_max -! read(iu_vmec,'(i4,9(a1,e15.7))') js_vmec(j),tb,r_vmec(j),tb, -! 1 roa_vmec(j),tb,chip_vmec(j),tb,psip_vmec(j),tb,btheta_vmec(j), -! 2 tb,bzeta_vmec(j),tb,vp_vmec(j),tb,bsq_vmec(j),tb,iota_vmec(j) -! end do -! close(iu_vmec) + top = 0.0 + bottom = 0.0 + + !Read VMEC wout data file + IF (.not. lwout_opened) CALL read_wout_file(TRIM(run_ident),ierr) ! Assign variables for the current surface chip = iota_vmec(js)*phipf_vmec(js); !note this only works for non-RFP - psip = phipf_vmec(js) - ! bsq = bsq_vmec(js) + psip = phipf_vmec(js) btheta = btheta_vmec(js); bzeta = bzeta_vmec(js) iota = iota_vmec(js); roa_surf = sqrt(phi_vmec(js)/phi_vmec(ns_vmec)) @@ -164,11 +185,12 @@ subroutine read_vmec_file_2(js,run_ident) vol_p = vp_vmec(js) Rmajor = Rmajor_vmec arad = Aminor_vmec - ! Now calc Bsq + + ! Now calc = int(B*B*sqrt(g)/sqrt(g)) bsq = 0.0 - DO v = 1, 360 + DO v = 1, nv_int zeta = TWOPI*REAL(v-1)/359. - DO u = 1, 360 + DO u = 1, nu_int theta = TWOPI*REAL(u-1)/359. bbf = 0.0; jacob_vmec = 0.0; DO mn = 1, mnmax_nyq @@ -182,10 +204,7 @@ subroutine read_vmec_file_2(js,run_ident) END DO END DO bsq = top/bottom - !deallocate variables -! deallocate(js_vmec, r_vmec, roa_vmec, chip_vmec, psip_vmec) -! deallocate(btheta_vmec, bzeta_vmec, vp_vmec, bsq_vmec,iota_vmec) - end subroutine read_vmec_file_2 + end subroutine read_vmec_file_wout c c---------------------------------------------------------------------------- c Reads the file "plasma_profiles*.dat" and assigns plasma parameters for the From 308d4ac5c500fcd4e45927bbddcbf8209eb3d238 Mon Sep 17 00:00:00 2001 From: lazersos Date: Wed, 28 Aug 2024 11:04:16 +0900 Subject: [PATCH 12/47] PENTA: Fixed bug --- PENTA/Sources/read_input_file_mod.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PENTA/Sources/read_input_file_mod.f b/PENTA/Sources/read_input_file_mod.f index bc204c093..c6394c73e 100644 --- a/PENTA/Sources/read_input_file_mod.f +++ b/PENTA/Sources/read_input_file_mod.f @@ -73,7 +73,7 @@ subroutine set_Utilde2_file(num_pts,roa_Utilde,U_tilde2, REAL(rknd), INTENT(OUT) :: U2 !dummy variables !local varibles - integer(iknd) :: j, js_surf, num_pts, kord_prof + integer(iknd) :: kord_prof real(rknd), dimension(:), allocatable :: U2_knot_array, spl_U2 !spline fit profile and evaluate at r/a of current surface From 536471ce8ff5e2c2fc115db1e4ff9d85512e0953 Mon Sep 17 00:00:00 2001 From: lazersos Date: Mon, 2 Sep 2024 14:04:32 +0900 Subject: [PATCH 13/47] PENTA: Added PENTA input module. --- PENTA/ObjectList | 1 + PENTA/PENTA.dep | 5 ++ PENTA/Sources/penta_functions_mod.f | 6 +- PENTA/Sources/penta_input_mod.f | 103 ++++++++++++++++++++++++++++ 4 files changed, 112 insertions(+), 3 deletions(-) create mode 100644 PENTA/Sources/penta_input_mod.f diff --git a/PENTA/ObjectList b/PENTA/ObjectList index 6a0b86914..7b008d112 100644 --- a/PENTA/ObjectList +++ b/PENTA/ObjectList @@ -1,4 +1,5 @@ ObjectFiles = \ +penta_input_mod.o \ penta_main.o \ penta_functions_mod.o \ penta_modules.o \ diff --git a/PENTA/PENTA.dep b/PENTA/PENTA.dep index 960e6f136..04320ae86 100644 --- a/PENTA/PENTA.dep +++ b/PENTA/PENTA.dep @@ -13,6 +13,11 @@ read_input_file_mod.o : \ bspline90_22.o \ penta_modules.o \ ../../LIBSTELL/$(LOCTYPE)/read_wout_mod.o + + +penta_input_mod.o : \ + penta_modules.o \ + ../../LIBSTELL/$(LOCTYPE)/safe_open_mod.o penta_subroutines.o : \ diff --git a/PENTA/Sources/penta_functions_mod.f b/PENTA/Sources/penta_functions_mod.f index 962270338..f94632a24 100644 --- a/PENTA/Sources/penta_functions_mod.f +++ b/PENTA/Sources/penta_functions_mod.f @@ -1,7 +1,7 @@ module penta_functions_mod - use penta_kind_mod - implicit none - contains + use penta_kind_mod + implicit none + contains cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc diff --git a/PENTA/Sources/penta_input_mod.f b/PENTA/Sources/penta_input_mod.f new file mode 100644 index 000000000..e4640084d --- /dev/null +++ b/PENTA/Sources/penta_input_mod.f @@ -0,0 +1,103 @@ + module penta_input_mod + use penta_kind_mod + implicit none + INTEGER(iknd), parameter :: num_ion_max = 20_iknd !Maximum number of ion species to be read from namelist file + INTEGER(iknd) :: num_ion_species ! Number of ion species + REAL(rknd), DIMENSION(num_ion_max) :: Z_ion_init,miomp_init + + namelist / ion_params / num_ion_species,Z_ion_init,miomp_init + contains + + SUBROUTINE init_ion_params_nml + IMPLICIT NONE + num_ion_species = 1 + Z_ion_init = 1.0_rknd + miomp_init = 1.0_rknd + END SUBROUTINE init_ion_params_nml + + SUBROUTINE read_ion_params_nml(filename,istat) + USE safe_open_mod + IMPLICIT NONE + ! Subroutine parameters + CHARACTER(256), INTENT(INOUT) :: filename + INTEGER(iknd), INTENT(INOUT) :: istat + ! Local Variables + LOGICAL :: lexist + INTEGER(iknd) :: iunit + CHARACTER(256) :: filename_local + + ! Handle filename + filename_local = TRIM(filename) + + !Check which file is available + INQUIRE(FILE=TRIM(filename_local),EXIST=lexist) + IF (.not. lexist) THEN + WRITE(6,*) ' -Could not open input file trying ion_params' + filename_local = 'ion_params' + INQUIRE(FILE=TRIM(filename_local),EXIST=lexist) + IF (.not. lexist) THEN + WRITE(6,*) ' -Could not open ion_params file' + istat = -1 + RETURN + ENDIF + ENDIF + + ! Open and read the file + CALL safe_open(iunit,istat,TRIM(filename_local),'old', + 1 'formatted') + IF (istat /= 0) THEN + WRITE(6,*) 'ERROR: Could no open: ',TRIM(filename_local) + RETURN + ENDIF + READ(iunit,NML=ion_params,IOSTAT=istat) + IF (istat /= 0) THEN + WRITE(6,*) 'ERROR: ion_params in ',TRIM(filename_local) + !backspace(iu_nl) + !read(iu_nl,fmt='(A)') line + !write(6,'(A)') 'Invalid line in namelist: '//TRIM(line) + RETURN + ENDIF + CLOSE(iunit) + RETURN + END SUBROUTINE read_ion_params_nml + + SUBROUTINE write_ion_params_nml(iunit) + IMPLICIT NONE + INTEGER(iknd), INTENT(inout) :: iunit + CHARACTER(LEN=*), PARAMETER :: outint = "(2X,A,1X,'=',1X,I0)" + INTEGER(iknd) :: k + WRITE(iunit,'(A)') '&ION_PARAMS' + WRITE(iunit,'(A)') + 1'----------------------------------------------------------------' + WRITE(iunit,outint) 'NUM_ION_SPECIES',num_ion_species + WRITE(iunit,"(2X,A,1X,'=',4(1X,ES22.12E3))") + 1'Z_ION_INIT',(Z_ion_init(k), k=1,num_ion_species) + WRITE(iunit,"(2X,A,1X,'=',4(1X,ES22.12E3))") + 1'MIOMP_INIT',(miomp_init(k), k=1,num_ion_species) + WRITE(iunit,'(A)') '/\n' + END SUBROUTINE write_ion_params_nml + + SUBROUTINE write_ion_params_namelist_byfile(filename) + USE safe_open_mod + IMPLICIT NONE + CHARACTER(LEN=*), INTENT(in) :: filename + INTEGER(iknd) :: iunit, istat + LOGICAL :: lexists + + iunit = 100 + istat = 0 + INQUIRE(FILE=TRIM(filename),exist=lexists) + IF (lexists) THEN + OPEN(unit=iunit, file=TRIM(filename), + 1 iostat=istat, status="old", position="append") + ELSE + OPEN(unit=iunit, file=TRIM(filename), + 1 iostat=istat, status="new") + END IF + IF (istat .ne. 0) RETURN + CALL write_ion_params_nml(iunit) + CLOSE(iunit) + + RETURN + END SUBROUTINE write_ion_params_namelist_byfile + end module penta_input_mod \ No newline at end of file From 9b7a5fb7b8c70bd8af9eeb266317ff356ad9e41c Mon Sep 17 00:00:00 2001 From: lazersos Date: Wed, 16 Oct 2024 13:48:57 +0200 Subject: [PATCH 14/47] LIBSTELL: Added thrift_globals.f90 for future python interfacing. --- LIBSTELL/LIBSTELL.dep | 3 + LIBSTELL/ObjectList | 1 + LIBSTELL/Sources/Modules/thrift_globals.f90 | 75 +++++++++++++++++++++ 3 files changed, 79 insertions(+) create mode 100644 LIBSTELL/Sources/Modules/thrift_globals.f90 diff --git a/LIBSTELL/LIBSTELL.dep b/LIBSTELL/LIBSTELL.dep index 5b09cc036..f7016a4b2 100644 --- a/LIBSTELL/LIBSTELL.dep +++ b/LIBSTELL/LIBSTELL.dep @@ -52,6 +52,9 @@ beams3d_input_mod.o: \ beams3d_globals.o : \ stel_kinds.o +thrift_globals.o : \ + stel_kinds.o + diagno_runtime.o : \ ezspline.o \ stel_kinds.o diff --git a/LIBSTELL/ObjectList b/LIBSTELL/ObjectList index a04acd792..79b36ba5c 100644 --- a/LIBSTELL/ObjectList +++ b/LIBSTELL/ObjectList @@ -1,4 +1,5 @@ ObjectFiles = \ +thrift_globals.o \ read_nescoil_mod.o \ stellopt_targets.o \ stellopt_vars.o \ diff --git a/LIBSTELL/Sources/Modules/thrift_globals.f90 b/LIBSTELL/Sources/Modules/thrift_globals.f90 new file mode 100644 index 000000000..5c252bc0f --- /dev/null +++ b/LIBSTELL/Sources/Modules/thrift_globals.f90 @@ -0,0 +1,75 @@ +!----------------------------------------------------------------------- +! Module: thrift_globals +! Authors: S. Lazerson (samuel.lazerson@gauss-fusion.com) +! Date: 10/16/2024 +! Description: This module contains the THRIFT global variables +! needed by the input namelist. +!----------------------------------------------------------------------- + MODULE thrift_globals +!----------------------------------------------------------------------- +! Libraries +!----------------------------------------------------------------------- + USE stel_kinds, ONLY: rprec + +!----------------------------------------------------------------------- +! Module Variables +!----------------------------------------------------------------------- + IMPLICIT NONE + + ! Moved from thrift_vars + LOGICAL :: lverbj + INTEGER :: nrho, ntimesteps, n_eq, npicard + REAL(rprec) :: tstart, tend, jtol, picard_factor, boot_factor + + ! Moved from thrift_vars (For ECCD in general) + INTEGER, PARAMETER :: ntime_ecrh = 200 + REAL(rprec), DIMENSION(ntime_ecrh) :: PECRH_AUX_T, PECRH_AUX_F + REAL(rprec) :: ecrh_rc, ecrh_w + + ! Moved from thrift_vars (for TRAVIS) + INTEGER, PARAMETER :: nsys = 16 + INTEGER :: nra_ecrh, nphi_ecrh + INTEGER, DIMENSION(nsys) :: wmode_ecrh + REAL(rprec), DIMENSION(nsys) :: freq_ecrh, power_ecrh + REAL(rprec), DIMENSION(nsys,3) :: antennaposition_ecrh, & + targetposition_ecrh, rbeam_ecrh, rfocus_ecrh + + ! Moved from thrift_vars (For DKES) + INTEGER, PARAMETER :: DKES_NS_MAX = 64 + INTEGER, PARAMETER :: DKES_NSTAR_MAX = 32 + INTEGER :: nruns_dkes + INTEGER, DIMENSION(:), POINTER :: DKES_rundex + INTEGER, DIMENSION(DKES_NS_MAX) :: DKES_K + REAL(rprec), DIMENSION(DKES_NSTAR_MAX) :: DKES_Erstar, DKES_Nustar + + ! Moved from thrift_runtime + INTEGER :: nparallel_runs, mboz, nboz + CHARACTER(256) :: bootstrap_type, eccd_type, vessel_ecrh, & + mirror_ecrh, targettype_ecrh, antennatype_ecrh + + + + CONTAINS + + ! These expose the global variables through ctypes + INTEGER FUNCTION getmaxtimeecrh() + IMPLICIT NONE + getmaxtimeecrh = ntime_ecrh + END FUNCTION getmaxtimeecrh + + INTEGER FUNCTION getmaxsys() + IMPLICIT NONE + getmaxsys = nsys + END FUNCTION getmaxsys + + INTEGER FUNCTION getmaxns() + IMPLICIT NONE + getmaxns = DKES_NS_MAX + END FUNCTION getmaxns + + INTEGER FUNCTION getmaxnstar() + IMPLICIT NONE + getmaxnstar =DKES_NSTAR_MAX + END FUNCTION getmaxnstar + + END MODULE thrift_globals From 81901ab1d98309ccb1bbf9ed89b3fc0c25d1cb1c Mon Sep 17 00:00:00 2001 From: lazersos Date: Wed, 16 Oct 2024 13:56:24 +0200 Subject: [PATCH 15/47] LIBSTELL: Added thrift_input_mod.f90 to LIBSTELL --- LIBSTELL/LIBSTELL.dep | 5 +++++ LIBSTELL/ObjectList | 1 + LIBSTELL/Sources/Modules/thrift_globals.f90 | 4 ++-- .../Sources/Modules}/thrift_input_mod.f90 | 14 ++++++++++---- 4 files changed, 18 insertions(+), 6 deletions(-) rename {THRIFT/Sources => LIBSTELL/Sources/Modules}/thrift_input_mod.f90 (95%) diff --git a/LIBSTELL/LIBSTELL.dep b/LIBSTELL/LIBSTELL.dep index f7016a4b2..ddd059785 100644 --- a/LIBSTELL/LIBSTELL.dep +++ b/LIBSTELL/LIBSTELL.dep @@ -55,6 +55,11 @@ beams3d_globals.o : \ thrift_globals.o : \ stel_kinds.o +thrift_input_mod.o: \ + stel_kinds.o \ + safe_open_mod.o \ + thrift_globals.o + diagno_runtime.o : \ ezspline.o \ stel_kinds.o diff --git a/LIBSTELL/ObjectList b/LIBSTELL/ObjectList index 79b36ba5c..829936646 100644 --- a/LIBSTELL/ObjectList +++ b/LIBSTELL/ObjectList @@ -1,4 +1,5 @@ ObjectFiles = \ +thrift_input_mod.o \ thrift_globals.o \ read_nescoil_mod.o \ stellopt_targets.o \ diff --git a/LIBSTELL/Sources/Modules/thrift_globals.f90 b/LIBSTELL/Sources/Modules/thrift_globals.f90 index 5c252bc0f..d60bbf5db 100644 --- a/LIBSTELL/Sources/Modules/thrift_globals.f90 +++ b/LIBSTELL/Sources/Modules/thrift_globals.f90 @@ -17,8 +17,8 @@ MODULE thrift_globals IMPLICIT NONE ! Moved from thrift_vars - LOGICAL :: lverbj - INTEGER :: nrho, ntimesteps, n_eq, npicard + LOGICAL :: lverbj, leccd, lnbcd, lohmic + INTEGER :: nrho, ntimesteps, n_eq, npicard, nsj REAL(rprec) :: tstart, tend, jtol, picard_factor, boot_factor ! Moved from thrift_vars (For ECCD in general) diff --git a/THRIFT/Sources/thrift_input_mod.f90 b/LIBSTELL/Sources/Modules/thrift_input_mod.f90 similarity index 95% rename from THRIFT/Sources/thrift_input_mod.f90 rename to LIBSTELL/Sources/Modules/thrift_input_mod.f90 index d52d86b2d..383214588 100644 --- a/THRIFT/Sources/thrift_input_mod.f90 +++ b/LIBSTELL/Sources/Modules/thrift_input_mod.f90 @@ -11,8 +11,7 @@ MODULE thrift_input_mod ! Libraries !----------------------------------------------------------------------- USE stel_kinds, ONLY: rprec - USE thrift_vars - USE thrift_runtime + USE thrift_globals USE safe_open_mod !----------------------------------------------------------------------- @@ -107,16 +106,23 @@ SUBROUTINE read_thrift_input(filename, istat) INQUIRE(FILE=TRIM(filename),EXIST=lexist) IF (.not.lexist) stop 'Could not find input file' CALL safe_open(iunit,istat,TRIM(filename),'old','formatted') - IF (istat /= 0) CALL handle_err(NAMELIST_READ_ERR,'thrift_input in: '//TRIM(filename),istat) + IF (istat /= 0) THEN + WRITE(6,'(A)') 'ERROR opening file: ',TRIM(filename) + CALL FLUSH(6) + STOP + END IF READ(iunit,NML=thrift_input,IOSTAT=istat) IF (istat /= 0) THEN + WRITE(6,'(A)') 'ERROR reading namelist THRIFT_INPUT from file: ',TRIM(filename) backspace(iunit) read(iunit,fmt='(A)') line write(6,'(A)') 'Invalid line in namelist: '//TRIM(line) - CALL handle_err(NAMELIST_READ_ERR,'thrift_input in: '//TRIM(filename),istat) + CALL FLUSH(6) + STOP END IF CLOSE(iunit) END IF + CALL tolower(bootstrap_type) CALL tolower(eccd_type) leccd = eccd_type .ne. '' From 41ff46d66222b8dfff3d3d7f66f089ff50070e3d Mon Sep 17 00:00:00 2001 From: lazersos Date: Wed, 16 Oct 2024 14:17:40 +0200 Subject: [PATCH 16/47] THRIFT: Cleanup after creating thrift_globals.f90 in LIBSTELL --- THRIFT/ObjectList | 1 - THRIFT/Sources/thrift_runtime.f90 | 25 +++++++++++----------- THRIFT/Sources/thrift_vars.f90 | 35 +++++++++---------------------- THRIFT/THRIFT.dep | 14 ++----------- 4 files changed, 25 insertions(+), 50 deletions(-) diff --git a/THRIFT/ObjectList b/THRIFT/ObjectList index 634831d20..1da248bf7 100644 --- a/THRIFT/ObjectList +++ b/THRIFT/ObjectList @@ -21,7 +21,6 @@ thrift_vars.o \ thrift_init.o \ thrift_init_mpisubgroup.o \ thrift_paraexe.o \ -thrift_input_mod.o \ thrift_profiles_mod.o \ thrift_interface_mod.o \ thrift_funcs.o \ diff --git a/THRIFT/Sources/thrift_runtime.f90 b/THRIFT/Sources/thrift_runtime.f90 index 67c737160..c6a4eab77 100644 --- a/THRIFT/Sources/thrift_runtime.f90 +++ b/THRIFT/Sources/thrift_runtime.f90 @@ -8,16 +8,20 @@ ! v0.00 11/XX/22 - Generally used to track major version information !----------------------------------------------------------------------- MODULE thrift_runtime - !----------------------------------------------------------------------- + !------------------------------------------------------------------- ! Libraries - !----------------------------------------------------------------------- + !------------------------------------------------------------------- USE stel_kinds, ONLY: rprec + USE thrift_globals, ONLY: nparallel_runs, mboz, nboz, & + bootstrap_type, eccd_type, vessel_ecrh, & + mirror_ecrh, targettype_ecrh, & + antennatype_ecrh USE mpi_params USE EZspline - !----------------------------------------------------------------------- + !------------------------------------------------------------------- ! Module Variables ! lverb Logical to control screen output - !---------------------------------------------------------------------- + !------------------------------------------------------------------- IMPLICIT NONE INTEGER, PARAMETER :: MPI_CHECK = 0 @@ -68,20 +72,17 @@ MODULE thrift_runtime DOUBLE PRECISION, PARAMETER :: e_charge = 1.60217662E-19 !e_c LOGICAL :: lverb, lvmec, lread_input, limas - INTEGER :: nprocs_thrift, nparallel_runs, mboz, nboz, ier_paraexe, & + INTEGER :: nprocs_thrift, ier_paraexe, & mytimestep, nsubsteps REAL(rprec) :: pi, pi2, invpi2, mu0, to3 - CHARACTER(256) :: id_string, prof_string, bootstrap_type, & - eccd_type, nbcd_type, & - proc_string, vessel_ecrh, mirror_ecrh, & - targettype_ecrh, antennatype_ecrh, & - magdiag_coil + CHARACTER(256) :: id_string, prof_string, & + nbcd_type, proc_string, magdiag_coil REAL(rprec), PARAMETER :: THRIFT_VERSION = 0.50 - !----------------------------------------------------------------------- + !------------------------------------------------------------------- ! Subroutines ! handle_err Controls Program Termination - !----------------------------------------------------------------------- + !------------------------------------------------------------------- CONTAINS SUBROUTINE handle_err(error_num, string_val, ierr) diff --git a/THRIFT/Sources/thrift_vars.f90 b/THRIFT/Sources/thrift_vars.f90 index 21290a09c..f78663517 100644 --- a/THRIFT/Sources/thrift_vars.f90 +++ b/THRIFT/Sources/thrift_vars.f90 @@ -10,6 +10,14 @@ MODULE thrift_vars ! Libraries !------------------------------------------------------------------- USE stel_kinds, ONLY: rprec + USE thrift_globals, ONLY: lverbj, nrho, ntimesteps, n_eq, npicard, & + tstart, tend, jtol, picard_factor, boot_factor, ntime_ecrh, & + pecrh_aux_t, pecrh_aux_f, ecrh_rc, ecrh_w, nsys, nra_ecrh, & + nphi_ecrh, wmode_ecrh, freq_ecrh, power_ecrh, & + antennaposition_ecrh, targetposition_ecrh, rbeam_ecrh, & + rfocus_ecrh, DKES_NS_MAX, DKES_NSTAR_MAX, nruns_dkes, & + DKES_rundex, DKES_K, dkes_Erstar, dkes_Nustar, nsj, leccd, & + lnbcd, lohmic !------------------------------------------------------------------- ! Module Variables ! leccd Calc Elec. Cyclo. Current Drive @@ -56,10 +64,9 @@ MODULE thrift_vars !------------------------------------------------------------------- IMPLICIT NONE - LOGICAL :: leccd, lnbcd, lohmic, ldiagno, lscreen_subcodes, lverbj + LOGICAL :: ldiagno, lscreen_subcodes LOGICAL, DIMENSION(:), ALLOCATABLE :: lbooz - INTEGER :: ntimesteps, nrho, nsj, npicard, n_eq,& - win_thrift_j,win_thrift_i,win_thrift_ugrid, & + INTEGER :: win_thrift_j,win_thrift_i,win_thrift_ugrid, & win_thrift_jplasma, win_thrift_iplasma, & win_thrift_jboot, win_thrift_iboot, & win_thrift_jeccd, win_thrift_ieccd, & @@ -77,7 +84,6 @@ MODULE thrift_vars win_thrift_alpha1, win_thrift_alpha2, win_thrift_alpha3, win_thrift_alpha4, & win_thrift_matld, win_thrift_matmd, win_thrift_matud, win_thrift_matrhs, & win_thrift_bvav - REAL(rprec) :: tstart, tend, jtol, picard_factor, boot_factor REAL(rprec), DIMENSION(:), POINTER :: THRIFT_RHO(:), THRIFT_RHOFULL(:), THRIFT_PHIEDGE(:), & THRIFT_S(:), THRIFT_SNOB(:), THRIFT_T(:) REAL(rprec), DIMENSION(:,:), POINTER :: & @@ -97,26 +103,5 @@ MODULE thrift_vars THRIFT_MATLD, THRIFT_MATMD, THRIFT_MATUD, THRIFT_MATRHS, & THRIFT_BVAV - ! For ECCD in general - INTEGER, PARAMETER :: ntime_ecrh = 200 - REAL(rprec), DIMENSION(ntime_ecrh) :: PECRH_AUX_T, PECRH_AUX_F - REAL(rprec) :: ecrh_rc, ecrh_w - - ! For TRAVIS - INTEGER, PARAMETER :: nsys = 16 - INTEGER :: nra_ecrh, nphi_ecrh - INTEGER, DIMENSION(nsys) :: wmode_ecrh - REAL(rprec), DIMENSION(nsys) :: freq_ecrh, power_ecrh - REAL(rprec), DIMENSION(nsys,3) :: antennaposition_ecrh, & - targetposition_ecrh,rbeam_ecrh,rfocus_ecrh - - ! For DKES - INTEGER, PARAMETER :: DKES_NS_MAX = 64 - INTEGER, PARAMETER :: DKES_NSTAR_MAX = 32 - INTEGER :: nruns_dkes - INTEGER, DIMENSION(:), POINTER :: DKES_rundex - INTEGER, DIMENSION(DKES_NS_MAX) :: DKES_K - REAL(rprec), DIMENSION(DKES_NSTAR_MAX) :: DKES_Erstar, DKES_Nustar - END MODULE thrift_vars diff --git a/THRIFT/THRIFT.dep b/THRIFT/THRIFT.dep index 38ac7e9a9..018053f37 100644 --- a/THRIFT/THRIFT.dep +++ b/THRIFT/THRIFT.dep @@ -1,7 +1,6 @@ thrift_bootsj.o : \ thrift_runtime.o \ thrift_vars.o \ - thrift_input_mod.o \ thrift_profiles_mod.o \ thrift_equil.o \ thrift_funcs.o \ @@ -27,7 +26,6 @@ thrift_boozer.o : \ thrift_dkes.o : \ thrift_runtime.o \ thrift_vars.o \ - thrift_input_mod.o \ thrift_profiles_mod.o \ thrift_equil.o \ thrift_funcs.o \ @@ -78,7 +76,6 @@ thrift_evolve.o : \ thrift_init.o : \ thrift_runtime.o \ - thrift_input_mod.o \ thrift_vars.o \ thrift_profiles_mod.o \ $(LIB_DIR)/$(LOCTYPE)/mpi_inc.o \ @@ -88,17 +85,10 @@ thrift_init.o : \ thrift_init_mpisubgroup.o : \ thrift_runtime.o \ - thrift_input_mod.o \ $(LIB_DIR)/$(LOCTYPE)/mpi_inc.o \ $(LIB_DIR)/$(LOCTYPE)/mpi_params.o \ $(LIB_DIR)/$(LOCTYPE)/safe_open_mod.o -thrift_input_mod.o : \ - thrift_runtime.o \ - thrift_vars.o \ - $(LIB_DIR)/$(LOCTYPE)/stel_kinds.o \ - $(LIB_DIR)/$(LOCTYPE)/safe_open_mod.o - thrift_interface_mod.o : \ thrift_runtime.o \ $(LIB_DIR)/$(LOCTYPE)/mpi_inc.o \ @@ -139,7 +129,6 @@ thrift_main.o : \ thrift_paraexe.o : \ thrift_runtime.o \ - thrift_input_mod.o \ $(LIB_DIR)/$(LOCTYPE)/mpi_inc.o \ $(LIB_DIR)/$(LOCTYPE)/mpi_params.o \ $(LIB_DIR)/$(LOCTYPE)/vmec_input.o \ @@ -166,6 +155,7 @@ thrift_reinit_vmec.o : \ thrift_runtime.o : \ $(LIB_DIR)/$(LOCTYPE)/stel_kinds.o \ + $(LIB_DIR)/$(LOCTYPE)/thrift_globals.o \ $(LIB_DIR)/$(LOCTYPE)/mpi_inc.o \ $(LIB_DIR)/$(LOCTYPE)/mpi_params.o @@ -200,7 +190,6 @@ thrift_run_nbcd.o : \ thrift_travis.o : \ thrift_runtime.o \ - thrift_input_mod.o \ thrift_vars.o \ thrift_equil.o \ thrift_profiles_mod.o \ @@ -210,6 +199,7 @@ thrift_travis.o : \ $(LIB_DIR)/$(LOCTYPE)/mpi_inc.o thrift_vars.o : \ + $(LIB_DIR)/$(LOCTYPE)/thrift_globals.o \ $(LIB_DIR)/$(LOCTYPE)/stel_kinds.o thrift_write.o : \ From 1871667c038bbe17a9999c7a5b0aa781d24d7291 Mon Sep 17 00:00:00 2001 From: lazersos Date: Wed, 16 Oct 2024 14:18:04 +0200 Subject: [PATCH 17/47] THRIFT: Added cleanup routines from STELLOPT. --- THRIFT/Sources/thrift_interface_mod.f90 | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/THRIFT/Sources/thrift_interface_mod.f90 b/THRIFT/Sources/thrift_interface_mod.f90 index 13416ede5..89d0e272a 100644 --- a/THRIFT/Sources/thrift_interface_mod.f90 +++ b/THRIFT/Sources/thrift_interface_mod.f90 @@ -77,8 +77,26 @@ END SUBROUTINE thrift_init_mpi_split SUBROUTINE thrift_cleanup IMPLICIT NONE + CHARACTER(200) :: cmdtxt = "" ! Clean up ierr_mpi = 0 + ! Remove files and cleanup directory +!DEC$ IF DEFINED (STELZIP) + IF (myid == master) THEN + ! Remove the *_opt* files + WRITE(6,*) ' Cleaning up files'; CALL FLUSH(6); ier = 0; ierr_mpi = 0; cmdtxt='' + CALL EXECUTE_COMMAND_LINE("rm -rf dcon* jxbout* mercier*",WAIT=.TRUE.,EXITSTAT=ier,CMDSTAT=ierr_mpi,CMDMSG=cmdtxt) + WRITE(6,*) ' rm: EXITSTAT=',ier,' CMDSTAT=',ierr_mpi; CALL FLUSH(6) + WRITE(6,*) ' MESSAGE: ',TRIM(cmdtxt); CALL FLUSH(6) + ! Zip up the results + WRITE(6,*) ' Zipping files'; CALL FLUSH(6); ier = 0; ierr_mpi = 0; cmdtxt='' + CALL EXECUTE_COMMAND_LINE("zip -r thrift_files.zip *",WAIT=.TRUE.,EXITSTAT=ier,CMDSTAT=ierr_mpi,CMDMSG=cmdtxt) + WRITE(6,*) ' zip: EXITSTAT=',ier,' CMDSTAT=',ierr_mpi; CALL FLUSH(6) + WRITE(6,*) ' MESSAGE: ',TRIM(cmdtxt); CALL FLUSH(6) + ier = 0; ierr_mpi=0 + END IF +!DEC$ ENDIF + !CALL thrift_free(MPI_COMM_SHARMEM) #if defined(MPI_OPT) !CALL MPI_BARRIER(MPI_COMM_MYWORLD,ierr_mpi) From 2bc3ba67e0969231fd82950a38ebddbc7e9416aa Mon Sep 17 00:00:00 2001 From: lazersos Date: Thu, 17 Oct 2024 10:29:04 +0200 Subject: [PATCH 18/47] PENTA: Added interface module to simplify interfacing to other codes. --- PENTA/ObjectList | 1 + PENTA/PENTA.dep | 12 + PENTA/Sources/penta_interface_mod.f90 | 877 ++++++++++++++++++++++++++ 3 files changed, 890 insertions(+) create mode 100644 PENTA/Sources/penta_interface_mod.f90 diff --git a/PENTA/ObjectList b/PENTA/ObjectList index ddd055cf3..c74604b38 100644 --- a/PENTA/ObjectList +++ b/PENTA/ObjectList @@ -1,4 +1,5 @@ ObjectFiles = \ +penta_interface_mod.o \ bspline.o \ coeff_var_pass.o \ io_unit_spec.o \ diff --git a/PENTA/PENTA.dep b/PENTA/PENTA.dep index 322013ea5..80b076008 100644 --- a/PENTA/PENTA.dep +++ b/PENTA/PENTA.dep @@ -20,6 +20,18 @@ penta.o : \ penta_math_routines_mod.o \ penta_subroutines.o +penta_interface_mod.o : \ + penta_kind_mod.o \ + io_unit_spec.o \ + read_input_file_mod.o \ + vmec_var_pass.o \ + pprof_pass.o \ + coeff_var_pass.o \ + phys_const.o \ + penta_functions_mod.o \ + penta_math_routines_mod.o \ + penta_subroutines.o + penta_functions_mod.o : \ penta_kind_mod.o \ phys_const.o \ diff --git a/PENTA/Sources/penta_interface_mod.f90 b/PENTA/Sources/penta_interface_mod.f90 new file mode 100644 index 000000000..ca80e2055 --- /dev/null +++ b/PENTA/Sources/penta_interface_mod.f90 @@ -0,0 +1,877 @@ +!----------------------------------------------------------------------- +! Module: PENTA_INTERFACE_MOD +! Authors: S. Lazerson +! Date: 10/16/2022 +! Description: Module provides a subroutine based interface to +! the PENTA code. +!----------------------------------------------------------------------- +MODULE PENTA_INTERFACE_MOD +!----------------------------------------------------------------------- +! Libraries +!----------------------------------------------------------------------- + USE penta_kind_mod + +!----------------------------------------------------------------------- +! Module Variables +!----------------------------------------------------------------------- + IMPLICIT NONE + + INTEGER(iknd), PARAMETER :: NUM_ROOTS_MAX = 10_iknd + INTEGER(iknd), PARAMETER :: NUM_ION_MAX = 20_iknd + + LOGICAL :: input_is_Er, log_interp, use_quanc8, read_U2_file, & + flux_cap, output_QoT_vs_Er, Add_Spitzer_to_D33, use_beam + INTEGER(iknd) :: num_Er_test, numKsteps, kord_pprof, keord, kcord, & + numargs, js, i_append, num_species, num_ion_species, Smax, & + iocheck, ie, ind_X, ind_A, ispec1, min_ind, iroot, num_roots + REAL(rknd) :: Kmin, Kmax, epsabs, epsrel, Er_min, Er_max, B_Eprl, & + U2, vth_e, loglambda, Er_test, abs_Er, min_Er, eaEr_o_kTe, cmin, & + cmax, emin, emax, sigma_par, sigma_par_Spitzer, J_BS + REAL(rknd), DIMENSION(NUM_ION_MAX) :: Z_ion_init, miomp_init + REAL(rknd), DIMENSION(NUM_ROOTS_MAX) :: Er_roots + REAL(rknd), DIMENSION(:), ALLOCATABLE :: ion_mass, Z_ion, vth_i, & + charges, dens, masses, temps, vths, dTdrs, dndrs, Xvec, Avec, & + Flows, Gammas, QoTs, xt_c, xt_e, Gamma_e_vs_Er, QoT_e_vs_Er, & + Er_test_vals, Jprl_ambi, sigma_par_ambi, sigma_par_Spitzer_ambi, & + J_BS_ambi + REAL(rknd), DIMENSION(:,:), ALLOCATABLE :: lmat, Dspl_D11, Dspl_D13, & + Dspl_D31, Dspl_D33, Dspl_Dex, Dspl_Dua, Dspl_Drat, Dspl_Drat2, & + Dspl_logD11, Dspl_logD33, cmesh, gamma_i_vs_er, QoT_i_vs_Er, & + Flows_ambi, gammas_ambi, QoTs_ambi, Jprl_parts, upol, utor + CHARACTER(LEN=10) :: Method + CHARACTER(LEN=100) :: arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, & + arg9, coeff_ext, run_ident, pprof_char, fpos, fstatus, str_num + +!----------------------------------------------------------------------- +! Module Namelists +!----------------------------------------------------------------------- + NAMELIST /ion_params/ num_ion_species, Z_ion_init, miomp_init + NAMELIST /run_params/ input_is_Er, log_interp, use_quanc8, & + read_U2_file, Add_Spitzer_to_D33, num_Er_test, numKsteps, & + kord_pprof, keord, kcord, Kmin, Kmax, epsabs, epsrel, Method, & + flux_cap, output_QoT_vs_Er, use_beam + +!----------------------------------------------------------------------- +! SUBROUTINES +!----------------------------------------------------------------------- + CONTAINS + + SUBROUTINE init_penta_input + IMPLICIT NONE + input_is_Er = .TRUE. + log_interp = .TRUE. + use_quanc8 = .FALSE. + read_U2_file = .TRUE. + flux_cap = .TRUE. + output_QoT_vs_Er = .FALSE. + Add_Spitzer_to_D33 = .TRUE. + use_beam = .FALSE. + num_Er_test = 50_iknd + numKsteps = 10000_iknd + kord_pprof = 3_iknd + keord = 2_iknd + kcord = 2_iknd + Kmin = 1.0E-5_rknd + Kmax = 2.0E+1_rknd + epsabs = 1.0E-8_rknd + epsrel = 1.0E-6_rknd + sigma_par = 0.0E+0_rknd + sigma_par_Spitzer = 0.0E+0_rknd + J_BS = 0.0E+0_rknd + method = 'DKES' + RETURN + END SUBROUTINE init_penta_input + + SUBROUTINE read_penta_ion_params_namelist(filename,istat) + USE safe_open_mod + IMPLICIT NONE + CHARACTER(*), INTENT(in) :: filename + INTEGER, INTENT(inout) :: istat + LOGICAL :: lexist + INTEGER :: iunit + CHARACTER(LEN=1000) :: line + istat = 0; iunit = 12 + INQUIRE(FILE=TRIM(filename),EXIST=lexist) + IF (.not.lexist) stop 'Could not find input file' + CALL safe_open(iunit,istat,TRIM(filename),'old','formatted') + IF (istat /= 0) THEN + WRITE(6,'(A)') 'ERROR opening file: ',TRIM(filename) + CALL FLUSH(6) + STOP + END IF + READ(iunit,NML=ion_params,IOSTAT=istat) + IF (istat /= 0) THEN + WRITE(6,'(A)') 'ERROR reading PENTA ion_params namelist from file: ',TRIM(filename) + backspace(iunit) + read(iunit,fmt='(A)') line + write(6,'(A)') 'Invalid line in namelist: '//TRIM(line) + CALL FLUSH(6) + STOP + END IF + CLOSE(iunit) + ! Electrons + ions + num_species = num_ion_species + 1_iknd + RETURN + END SUBROUTINE read_penta_ion_params_namelist + + SUBROUTINE read_penta_run_params_namelist(filename,istat) + USE safe_open_mod + IMPLICIT NONE + CHARACTER(*), INTENT(in) :: filename + INTEGER, INTENT(inout) :: istat + LOGICAL :: lexist + INTEGER :: iunit + CHARACTER(LEN=1000) :: line + istat = 0; iunit = 12 + INQUIRE(FILE=TRIM(filename),EXIST=lexist) + IF (.not.lexist) stop 'Could not find input file' + CALL safe_open(iunit,istat,TRIM(filename),'old','formatted') + IF (istat /= 0) THEN + WRITE(6,'(A)') 'ERROR opening file: ',TRIM(filename) + CALL FLUSH(6) + STOP + END IF + READ(iunit,NML=run_params,IOSTAT=istat) + IF (istat /= 0) THEN + WRITE(6,'(A)') 'ERROR reading PENTA ion_params namelist from file: ',TRIM(filename) + backspace(iunit) + read(iunit,fmt='(A)') line + write(6,'(A)') 'Invalid line in namelist: '//TRIM(line) + CALL FLUSH(6) + STOP + END IF + CLOSE(iunit) + RETURN + END SUBROUTINE read_penta_run_params_namelist + + SUBROUTINE penta_init_commandline + IMPLICIT NONE + CALL GETCARG(1, arg1, numargs) + IF (numargs .NE. 9) THEN + Write(6,'(A)') 'ERROR Incorrect number of input arguments, see penta.f90 for details' + CALL FLUSH(6) + STOP + END IF + CALL GETCARG(2, arg2, numargs) + CALL GETCARG(3, arg3, numargs) + CALL GETCARG(4, arg4, numargs) + CALL GETCARG(5, arg5, numargs) + CALL GETCARG(6, arg6, numargs) + CALL GETCARG(7, arg7, numargs) + CALL GETCARG(8, arg8, numargs) + CALL GETCARG(9, arg9, numargs) + + READ(arg2,*) Er_min + READ(arg3,*) Er_max + READ(arg4,*) js + READ(arg5,*) i_append + READ(arg8,*) B_Eprl + READ(arg9,*) Smax + coeff_ext = TRIM(ADJUSTL(arg1)) + run_ident = TRIM(ADJUSTL(arg6)) + pprof_char = TRIM(ADJUSTL(arg7)) + RETURN + END SUBROUTINE penta_init_commandline + + SUBROUTINE penta_allocate_species + USE pprof_pass, ONLY: ni, Ti, dnidr, dTidr + IMPLICIT NONE + ALLOCATE(ni(num_ion_species)) + ALLOCATE(Ti(num_ion_species)) ! Ion profile info + ALLOCATE(dnidr(num_ion_species)) + ALLOCATE(dTidr(num_ion_species)) + ALLOCATE(Z_ion(num_ion_species)) + ALLOCATE(ion_mass(num_ion_species)) ! Ion parameters + ALLOCATE(vth_i(num_ion_species)) + ALLOCATE(charges(num_species)) ! Parameters for + ALLOCATE(dens(num_species)) ! all species + ALLOCATE(vths(num_species)) + ALLOCATE(masses(num_species)) + ALLOCATE(Temps(num_species)) + ALLOCATE(dTdrs(num_species)) + ALLOCATE(dndrs(num_species)) + ALLOCATE(lmat((Smax+1)*num_species,(Smax+1)*num_species)) ! Clas. fric.coeffs + ALLOCATE(Xvec(num_species*2+1)) + ALLOCATE(Avec(num_species*3)) ! Thermo. Force vecs. + ALLOCATE(Flows((Smax+1)*num_species)) ! Prl flow moments + ALLOCATE(Gammas(num_species)) ! Rad fluxes + ALLOCATE(QoTs(num_species)) ! Rad energy fluxes + ALLOCATE(Gamma_i_vs_Er(num_Er_test,num_ion_species)) ! Ion flux vs Er + ALLOCATE(Gamma_e_vs_Er(num_Er_test)) ! Electron flux vs Er + ALLOCATE(Er_test_vals(num_Er_test)) ! Er to loop over + IF ( output_QoT_vs_Er ) THEN + ALLOCATE(QoT_i_vs_Er(num_Er_test,num_ion_species)) ! Ion flux vs Er + ALLOCATE(QoT_e_vs_Er(num_Er_test)) ! Electron flux vs Er + ENDIF + RETURN + END SUBROUTINE penta_allocate_species + + SUBROUTINE penta_allocate_dkescoeff + USE coeff_var_pass, ONLY: num_c, num_e + IMPLICIT NONE + ALLOCATE(xt_c(num_c + kcord)) + ALLOCATE(xt_e(num_e + keord)) + ALLOCATE(Dspl_D11(num_c,num_e)) + ALLOCATE(Dspl_D13(num_c,num_e)) + ALLOCATE(Dspl_D31(num_c,num_e)) + ALLOCATE(Dspl_D33(num_c,num_e)) + ALLOCATE(Dspl_Dex(num_c,num_e)) + ALLOCATE(Dspl_Drat(num_c,num_e)) + ALLOCATE(Dspl_Drat2(num_c,num_e)) + ALLOCATE(Dspl_DUa(num_c,num_e)) + ALLOCATE(Dspl_logD11(num_c,num_e)) + ALLOCATE(Dspl_logD33(num_c,num_e)) + ALLOCATE(cmesh(num_c,num_e)) + RETURN + END SUBROUTINE penta_allocate_dkescoeff + + SUBROUTINE penta_deallocate_species + USE pprof_pass, ONLY: ni, Ti, dnidr, dTidr + IMPLICIT NONE + IF (ALLOCATED(ni)) DEALLOCATE(ni) + IF (ALLOCATED(ti)) DEALLOCATE(ti) + IF (ALLOCATED(dnidr)) DEALLOCATE(dnidr) + IF (ALLOCATED(dTidr)) DEALLOCATE(dTidr) + IF (ALLOCATED(Z_ion)) DEALLOCATE(Z_ion) + IF (ALLOCATED(ion_mass)) DEALLOCATE(ion_mass) + IF (ALLOCATED(vth_i)) DEALLOCATE(vth_i) + IF (ALLOCATED(charges)) DEALLOCATE(charges) + IF (ALLOCATED(dens)) DEALLOCATE(dens) + IF (ALLOCATED(vths)) DEALLOCATE(vths) + IF (ALLOCATED(masses)) DEALLOCATE(masses) + IF (ALLOCATED(Temps)) DEALLOCATE(Temps) + IF (ALLOCATED(dTdrs)) DEALLOCATE(dTdrs) + IF (ALLOCATED(dndrs)) DEALLOCATE(dndrs) + IF (ALLOCATED(lmat)) DEALLOCATE(lmat) + IF (ALLOCATED(Xvec)) DEALLOCATE(Xvec) + IF (ALLOCATED(Avec)) DEALLOCATE(Avec) + IF (ALLOCATED(Flows)) DEALLOCATE(Flows) + IF (ALLOCATED(Gammas)) DEALLOCATE(Gammas) + IF (ALLOCATED(QoTs)) DEALLOCATE(QoTs) + IF (ALLOCATED(Gamma_i_vs_Er)) DEALLOCATE(Gamma_i_vs_Er) + IF (ALLOCATED(Gamma_e_vs_Er)) DEALLOCATE(Gamma_e_vs_Er) + IF (ALLOCATED(Er_test_vals)) DEALLOCATE(Er_test_vals) + IF (ALLOCATED(QoT_i_vs_Er)) DEALLOCATE(QoT_i_vs_Er) + IF (ALLOCATED(QoT_e_vs_Er)) DEALLOCATE(QoT_e_vs_Er) + RETURN + END SUBROUTINE penta_deallocate_species + + SUBROUTINE penta_deallocate_dkescoeff + IMPLICIT NONE + IF (ALLOCATED(xt_c)) DEALLOCATE(xt_c) + IF (ALLOCATED(xt_e)) DEALLOCATE(xt_e) + IF (ALLOCATED(Dspl_D11)) DEALLOCATE(Dspl_D11) + IF (ALLOCATED(Dspl_D13)) DEALLOCATE(Dspl_D13) + IF (ALLOCATED(Dspl_D31)) DEALLOCATE(Dspl_D31) + IF (ALLOCATED(Dspl_D33)) DEALLOCATE(Dspl_D33) + IF (ALLOCATED(Dspl_Dex)) DEALLOCATE(Dspl_Dex) + IF (ALLOCATED(Dspl_Drat)) DEALLOCATE(Dspl_Drat) + IF (ALLOCATED(Dspl_Drat2)) DEALLOCATE(Dspl_Drat2) + IF (ALLOCATED(Dspl_DUa)) DEALLOCATE(Dspl_DUa) + IF (ALLOCATED(Dspl_logD11)) DEALLOCATE(Dspl_logD11) + IF (ALLOCATED(Dspl_logD33)) DEALLOCATE(Dspl_logD33) + IF (ALLOCATED(cmesh)) DEALLOCATE(cmesh) + RETURN + END SUBROUTINE penta_deallocate_dkescoeff + + SUBROUTINE penta_read_input_files + USE vmec_var_pass, ONLY: roa_surf, arad, Bsq + USE pprof_pass + USE coeff_var_pass, ONLY: D11_mat, cmul, num_c + USE phys_const, ONLY: p_mass, elem_charge, e_mass + USE read_input_file_mod + USE PENTA_subroutines, ONLY : define_friction_coeffs + IMPLICIT NONE + CALL read_vmec_file_2(js,run_ident) + CALL read_pprof_file(pprof_char,num_ion_species,roa_surf,arad,kord_pprof) + CALL read_dkes_star_files(coeff_ext,Add_Spitzer_to_D33,Bsq) + IF (use_beam) THEN + CALL read_beam_file(roa_surf,kord_pprof) + ELSE + beam_force = 0._rknd + ENDIF + ! Optionally read file containing info. Else this is + ! calculated from the D11* coefficient at high nu/v and Er=0. + IF ( read_U2_file ) THEN + CALL read_Utilde2_file(roa_surf,U2,kord_pprof) + ELSE + U2=1.5d0*D11_mat(num_c,1)/cmul(num_c); + ENDIF + ! Change Er test range to V/cm if necessary + IF ( input_is_Er) THEN + Er_min = Er_min * Te / arad + Er_max = Er_max * Te / arad + ENDIF + ! Assign ion parameters + Z_ion = Z_ion_init(1:num_ion_species) + ion_mass = miomp_init(1:num_ion_species) * p_mass + ! Calculate thermal velocities + vth_i = Dsqrt(2._rknd*Ti*elem_charge/ion_mass) + vth_e = Dsqrt(2._rknd*Te*elem_charge/e_mass) + ! Calculate Coulomb logarithm + IF ( Te > 50._rknd ) THEN + loglambda = 25.3_rknd - 1.15_rknd*Dlog10(ne/1.e6_rknd) + 2.3_rknd*Dlog10(Te) + ELSE + loglambda = 23.4_rknd - 1.15_rknd*Dlog10(ne/1.e6_rknd) + 3.45_rknd*Dlog10(Te) + ENDIF + ! Assign arrays of parameters for all species (charge, mass, n, T, v_th, dTdr) + charges=elem_charge*(/-1._rknd,Z_ion/) + dens=(/ne, ni/) + masses=(/e_mass, ion_mass/) + Temps=(/Te,Ti/) + vths=(/vth_e,vth_i/) + dTdrs=(/dTedr,dTidr/) + dndrs=(/dnedr,dnidr/) + ! Define matrix of friction coefficients (lmat) + Call define_friction_coeffs(masses,charges,vths,Temps,dens,loglambda, & + num_species,Smax,lmat) + RETURN + END SUBROUTINE penta_read_input_files + + SUBROUTINE penta_fit_DXX_coef + USE coeff_var_pass + USE PENTA_subroutines, ONLY : fit_coeffs + IMPLICIT NONE + + ! Calculate fitting parameters to the D##* coefficients + Call fit_coeffs(cmul,efield,num_c,num_e,D11_mat,log_interp, & + kcord,keord,xt_c,xt_e,Dspl_D11,cmin,cmax,emin,emax) + Call fit_coeffs(cmul,efield,num_c,num_e,D13_mat,log_interp, & + kcord,keord,xt_c,xt_e,Dspl_D13,cmin,cmax,emin,emax) + Call fit_coeffs(cmul,efield,num_c,num_e,D31_mat,log_interp, & + kcord,keord,xt_c,xt_e,Dspl_D31,cmin,cmax,emin,emax) + Call fit_coeffs(cmul,efield,num_c,num_e,D33_mat,log_interp, & + kcord,keord,xt_c,xt_e,Dspl_D33,cmin,cmax,emin,emax) + + ! Fit log(D*) for D11 and D33 + Call fit_coeffs(cmul,efield,num_c,num_e,LOG(D11_mat),log_interp, & + kcord,keord,xt_c,xt_e,Dspl_logD11,cmin,cmax,emin,emax) + Call fit_coeffs(cmul,efield,num_c,num_e,LOG(D33_mat),log_interp, & + kcord,keord,xt_c,xt_e,Dspl_logD33,cmin,cmax,emin,emax) + RETURN + END SUBROUTINE penta_fit_DXX_coef + + SUBROUTINE penta_screen_info + IMPLICIT NONE + If ( i_append == 0 ) Then + WRITE(6,'(A)') "" + WRITE(6,'(A)') "Welcome to PENTA3, please note the following settings:" + WRITE(6,'(A)') + WRITE(6,'(A,I3)') ' Number of ion species: ',num_ion_species + IF ( input_is_Er ) THEN + WRITE(6,'(A)') 'Interpreting input range as Er (V/cm)' + ELSE + WRITE(6,'(A)') 'Interpreting input range as eEr/kTe' + ENDIF + IF ( log_interp ) THEN + WRITE(6,'(A)') 'Performing logarithmic interpolation in Er, cmul' + ELSE + WRITE(6,'(A)') 'Performing linear interpolation in Er,cmul' + ENDIF + IF ( use_quanc8 ) THEN + WRITE(6,'(A,2(A,E10.4))') & + ' Using quanc8 integrator with tolerances: ', & + 'abs: ',epsabs,' rel: ', epsrel + ELSE + WRITE(6,'(A,i6,A)') ' Using ',numKsteps, & + ' point integral approximation' + ENDIF + WRITE(6,'(a,2(" ",e15.4))') ' K range on convolution integral: ', & + Kmin, Kmax + IF ( Add_Spitzer_to_D33) & + WRITE(6,'(A)') 'Adding collisional (Spitzer) portion to D33 coefficient' + IF ( flux_cap ) & + WRITE(6,'(A)') 'Enforcing minimum radial diffusion coefficient = 0' + WRITE(6,'(A,I2)') ' Number of terms in Sonine expansion: ', Smax+1 + SELECT CASE (Method) + CASE ('T') + WRITE(6,'(A)') 'Using Taguchi Method' + CASE ('SN') + WRITE(6,'(A)') 'Using Sugama-Nishimura Method' + CASE ('MBT') + WRITE(6,'(A)') 'Using Maassberg-Beidler-Turkin Method' + CASE ('DKES') + WRITE(6,'(A)') 'Using DKES Method' + CASE DEFAULT + Write(6,'(3A)') ' Error: ''', Trim(Adjustl(Method)), & + ''' is not a valid Method' + STOP 'Error: Exiting, method select error in penta.f90 (1)' + END SELECT + WRITE(6,'(A)') "" + WRITE(6,'(A)') " /"," Er root(s) (V/cm)" + ENDIF + RETURN + END SUBROUTINE penta_screen_info + + SUBROUTINE penta_open_output + USE safe_open_mod + USE io_unit_spec + IMPLICIT NONE + INTEGER :: istat + istat = 0 + ! Set write status + IF (i_append == 0) THEN + fstatus = "unknown" + fpos = "asis" + ELSEIF (i_append == 1) THEN + fstatus = "old" + fpos = "append" + ELSE + Write(6,'(A)') 'PENTA: Bad value for i_append (0 or 1 expected)' + CALL FLUSH(6) + STOP 'Error: Exiting, i_append error in penta.f90' + END IF + ! Open files + !Open(unit=iu_flux_out, file="fluxes_vs_roa", position=Trim(Adjustl(fpos)),status=Trim(Adjustl(fstatus))) + CALL safe_open(iu_flux_out, istat, "fluxes_vs_roa", & + Trim(Adjustl(fstatus)), 'formatted',& + access_in=Trim(Adjustl(fpos))) + CALL safe_open(iu_pprof_out, istat, "plasma_profiles_check", & + Trim(Adjustl(fstatus)), 'formatted',& + access_in=Trim(Adjustl(fpos))) + CALL safe_open(iu_fvEr_out, istat, "fluxes_vs_Er", & + Trim(Adjustl(fstatus)), 'formatted',& + access_in=Trim(Adjustl(fpos))) + CALL safe_open(iu_flows_out, istat, "flows_vs_roa", & + Trim(Adjustl(fstatus)), 'formatted',& + access_in=Trim(Adjustl(fpos))) + CALL safe_open(iu_flowvEr_out, istat, "flows_vs_Er", & + Trim(Adjustl(fstatus)), 'formatted',& + access_in=Trim(Adjustl(fpos))) + CALL safe_open(iu_Jprl_out, istat, "Jprl_vs_roa", & + Trim(Adjustl(fstatus)), 'formatted',& + access_in=Trim(Adjustl(fpos))) + CALL safe_open(iu_contraflows_out, istat, "ucontra_vs_roa", & + Trim(Adjustl(fstatus)), 'formatted',& + access_in=Trim(Adjustl(fpos))) + IF (method == 'SN') & + CALL safe_open(iu_sigmas_out, istat, "sigmas_vs_roa", & + Trim(Adjustl(fstatus)), 'formatted',& + access_in=Trim(Adjustl(fpos))) + IF (output_QoT_vs_Er) THEN + CALL safe_open(iu_QoTvEr_out, istat, "QoTs_vs_Er", & + Trim(Adjustl(fstatus)), 'formatted',& + access_in=Trim(Adjustl(fpos))) + Write(iu_QoTvEr_out,'("*",/,"r/a Er[V/cm] Q_e/T_e [m**-2s**-1] ",& + " Q_i/T_i [m**-2s**-1]")') + END IF + ! WRITE Headers + IF (i_append == 0) THEN + ! Fluxes vs r/a + Write(iu_flux_out,'("*",/,"r/a Er[V/cm] eEr/kTe ", & + "Gamma_e [m**-2s**-1] Q_e/T_e [m**-2s**-1] ", & + "Gamma_i [m**-2s**-1] Q_i/T_i [m**-2s**-1]")') + ! Flows vs r/a + Write(iu_flows_out,'("*",/,"r/a Er[V/cm] eEr/kTe ", & + " / [m/sT] / [m/sT]")') + ! Plasma profile check + Write(iu_pprof_out,'("*",/,"r/a Te [eV] ne [m**-3] ", & + "dnedr [m**-4] dTedr [eV/m] Ti [eV] ni [m**-3] ", & + "dnidr [m**-4] dTidr [eV/m]")') + Write(iu_Jprl_out,'("*",/,"r/a Er [V/cm] eEr/kTe ", & + "Jprl_e [A/m**2] Jprli [A/m**2] Jprl [A/m**2] J_BS [A/m**2]")') + Write(iu_contraflows_out,'("*",/,"r/a Er [V/cm] eEr/kTe ", & + " [1/s] [1/s] ", & + " [1/s] [1/s]")') + ! Sigmas vs r/a + IF (Method == 'SN') & + Write(iu_sigmas_out,'("*",/,"r/a Er[V/cm] sigma_par [1/Ohm.m] ", & + " sigma_par_Spitzer [1/Ohm.m]")') + ! Legend for fluxes vs Er + Write(iu_fvEr_out,'("*",/,"r/a Er[V/cm] Gamma_e [m**-2s**-1] ",& + " Gamma_i [m**-2s**-1]")') + ! Legend for flows vs Er + Write(iu_flowvEr_out,'("*",/,"r/a Er[V/cm] ", & + " / [m/sT] / [m/sT]")') + END IF + RETURN + END SUBROUTINE penta_open_output + + SUBROUTINE penta_fit_rad_trans + USE coeff_var_pass + USE vmec_var_pass + USE PENTA_subroutines, ONLY : fit_coeffs + IMPLICIT NONE + ! Fit radial transport coefficients specific to different methods + SELECT CASE (Method) + CASE ('T', 'MBT') + cmesh = Spread(cmul,2,num_e) + ! Calculate the D11 coefficient minus the P-S contribution (Dex) + ! Also, do not allow for negative coefficients + CALL fit_coeffs(cmul,efield,num_c,num_e, & + Max(D11_mat-(2._rknd/3._rknd)*cmesh*U2,0._rknd), & + log_interp,kcord,keord,xt_c,xt_e,Dspl_Dex, & + cmin,cmax,emin,emax) + CASE ('SN') + ! Calculate fits to D31*/D33* (Drat) + CALL fit_coeffs(cmul,efield,num_c,num_e, & + D31_mat/D33_mat, & + log_interp,kcord,keord,xt_c,xt_e,Dspl_Drat, & + cmin,cmax,emin,emax) + ! Calculate fits to (D31*)**2/D33* (Drat2) + CALL fit_coeffs(cmul,efield,num_c,num_e, & + D31_mat*D31_mat/D33_mat, & + log_interp,kcord,keord,xt_c,xt_e,Dspl_Drat2, & + cmin,cmax,emin,emax) + cmesh = Spread(cmul,2,num_e) + ! Calculate coefficient for Ua term (DUa) + CALL fit_coeffs(cmul,efield,num_c,num_e, & + (2._rknd*Bsq/(3._rknd*D33_mat) - cmesh), & + log_interp,kcord,keord,xt_c,xt_e,Dspl_DUa, & + cmin,cmax,emin,emax) + ! Calculate coefficient for radial flux (Capped term) (Dex) + ! Also, do not allow for negative coefficients + CALL fit_coeffs(cmul,efield,num_c,num_e, & + Max(D11_mat-(2._rknd/3._rknd)*cmesh*U2+D31_mat*D31_mat/D33_mat, & + 0._rknd),log_interp,kcord,keord,xt_c,xt_e,Dspl_Dex, & + cmin,cmax,emin,emax) + CASE ('DKES') + CASE DEFAULT + WRITE(6,'(3a)') ' Error: ''', Trim(Adjustl(Method)), & + ''' is not a valid Method' + STOP 'Error: Exiting, method select error in penta.f90 (2)' + ENDSELECT + RETURN + END SUBROUTINE penta_fit_rad_trans + + SUBROUTINE penta_run_1_init + IMPLICIT NONE + INTEGER :: istat + istat = 0 + CALL read_penta_ion_params_namelist('ion_params',istat) + istat = 0 + CALL read_penta_run_params_namelist('run_params',istat) + CALL penta_init_commandline + CALL penta_allocate_species + CALL penta_read_input_files + CALL penta_screen_info + CALL penta_allocate_dkescoeff + CALL penta_fit_DXX_coef + CALL penta_open_output + CALL penta_fit_rad_trans + RETURN + END SUBROUTINE penta_run_1_init + + SUBROUTINE penta_run_2_efield + USE vmec_var_pass + USE pprof_pass + USE phys_const + USE io_unit_spec + USE coeff_var_pass + USE penta_math_routines_mod, ONLY: rlinspace + USE penta_functions_mod + USE PENTA_subroutines, ONLY: form_xvec + IMPLICIT NONE + ! Define array of Er values to test [V/m] + Er_test_vals = rlinspace(Er_min,Er_max,num_Er_test)*100._rknd + + ! Check for Er=0, doesn't work for log interpolation + min_Er = Minval(Dabs(Er_test_vals),DIM=1) + If ((log_interp .EQV. .true. ) .AND. ( Dabs(min_Er) <= elem_charge )) Then + min_ind = Minloc(Dabs(Er_test_vals),DIM=1) + If ( min_ind == Num_Er_test ) Then + Er_test_vals(min_ind) = Er_test_vals(min_ind - 1)/2._rknd + Else + Er_test_vals(min_ind) = Er_test_vals(min_ind + 1)/2._rknd + EndIf + Write(*,'(a,i4,a,f10.3)') 'Cannot use Er=0 with log_interp, using Er(', & + min_ind, ') = ', Er_test_vals(min_ind) + EndIf + ! Loop over Er to get fluxes as a function of Er + Do ie = 1,num_Er_test + Er_test = Er_test_vals(ie) + abs_Er = Abs(Er_test) + + ! Form thermodynamic force vector (Xvec) + Call form_Xvec(Er_test,Z_ion,B_Eprl,num_ion_species,Xvec) + + ! Form alternate thermodynamic force vector (Avec) + Do ispec1 = 1,num_species + ind_X = (ispec1-1)*2 + 1 + ind_A = (ispec1-1)*3 + 1 + + Avec(ind_A) = -Xvec(ind_X) / (Temps(ispec1)*elem_charge) & + - 2.5_rknd*dTdrs(ispec1)/Temps(ispec1) + Avec(ind_A+1) = -Xvec(ind_X+1) / (Temps(ispec1)*elem_charge) + Avec(ind_A+2) = Xvec(num_species*2+1)*charges(ispec1) & + * B0/(Temps(ispec1)*elem_charge*Sqrt(Bsq)) + & + beam_force/(Temps(ispec1)*elem_charge*dens(ispec1)) + Enddo + + ! Select the appropriate algorithm and calculate the flows and fluxes + SELECT CASE (Method) + Case ('T', 'MBT') + ! Calculate array of parallel flow moments + Flows = calc_flows_T(num_species,Smax,abs_Er,Temps,dens,vths,charges, & + masses,loglambda,B0,use_quanc8,Kmin,Kmax,numKsteps,log_interp,cmin, & + cmax,emin,emax,xt_c,xt_e,Dspl_D31,Dspl_logD33,num_c,num_e,kcord, & + keord,Avec,Bsq,lmat,J_BS) + Gammas = calc_fluxes_MBT(num_species,Smax,abs_Er,Temps,dens,vths, & + charges,masses,dTdrs,dndrs,loglambda,use_quanc8,Kmin,Kmax,numKsteps, & + log_interp,cmin,cmax,emin,emax,xt_c,xt_e,Dspl_logD11,Dspl_D31, & + Dspl_Dex,num_c,num_e,kcord,keord,Avec,lmat,Flows,U2,B0,flux_cap) + If ( output_QoT_vs_Er .EQV. .true. ) Then + QoTs = calc_QoTs_MBT(num_species,Smax,abs_Er,Temps,dens,vths,charges, & + masses,dTdrs,dndrs,loglambda,use_quanc8,Kmin,Kmax,numKsteps, & + log_interp,cmin,cmax,emin,emax,xt_c,xt_e,Dspl_logD11,Dspl_D31, & + Dspl_Dex,num_c,num_e,kcord,keord,Avec,lmat,Flows,U2,B0,flux_cap) + Endif + Case ('SN') + Flows = calc_flows_SN(num_species,Smax,abs_Er,Temps,dens,vths,charges, & + masses,loglambda,B0,use_quanc8,Kmin,Kmax,numKsteps,log_interp, & + cmin,cmax,emin,emax,xt_c,xt_e,Dspl_Drat,Dspl_DUa,num_c,num_e,kcord, & + keord,Avec,lmat,sigma_par,sigma_par_Spitzer,J_BS) + Gammas = calc_fluxes_SN(num_species,Smax,abs_Er,Temps,dens,vths,charges,& + masses,loglambda,use_quanc8,Kmin,Kmax,numKsteps,log_interp,cmin,cmax, & + emin,emax,xt_c,xt_e,Dspl_Drat,Dspl_Drat2,Dspl_Dex,Dspl_logD11, & + Dspl_D31,num_c,num_e,kcord,keord,Avec,Bsq,lmat,Flows,U2,dTdrs, & + dndrs,flux_cap) + If ( output_QoT_vs_Er .EQV. .true. ) Then + QoTs = calc_QoTs_SN(num_species,Smax,abs_Er,Temps,dens,vths,charges, & + masses,loglambda,use_quanc8,Kmin,Kmax,numKsteps,log_interp,cmin, & + cmax,emin,emax,xt_c,xt_e,Dspl_Drat,Dspl_Drat2,Dspl_Dex,Dspl_logD11, & + Dspl_D31,num_c,num_e,kcord,keord,Avec,Bsq,lmat,Flows,U2,dTdrs, & + dndrs,flux_cap) + Endif + Case ('DKES') + Flows = calc_flows_DKES(num_species,Smax,abs_Er,Temps,dens,vths,charges,& + masses,loglambda,B0,use_quanc8,Kmin,Kmax,numKsteps,log_interp,cmin, & + cmax,emin,emax,xt_c,xt_e,Dspl_D31,Dspl_logD33,num_c,num_e,kcord, & + keord,Avec,J_BS) + Gammas = calc_fluxes_DKES(num_species,abs_Er,Temps,dens,vths,charges, & + masses,loglambda,use_quanc8,Kmin,Kmax,numKsteps,log_interp,cmin,cmax, & + emin,emax,xt_c,xt_e,Dspl_logD11,Dspl_D31,num_c,num_e,kcord,keord, & + Avec,B0) + If ( output_QoT_vs_Er .EQV. .true. ) Then + QoTs = calc_QoTs_DKES(num_species,abs_Er,Temps,dens,vths,charges, & + masses,loglambda,use_quanc8,Kmin,Kmax,numKsteps,log_interp,cmin, & + cmax,emin,emax,xt_c,xt_e,Dspl_logD11,Dspl_D31,num_c,num_e,kcord, & + keord,Avec,B0) + Endif + Case Default + Write(6,'(3a)') ' Error: ''', Trim(Adjustl(Method)), & + ''' is not a valid Method' + Stop 'Error: Exiting, method select error in penta.f90 (3)' + END SELECT + + Gamma_e_vs_Er(ie) = Gammas(1) + Gamma_i_vs_Er(ie,:) = Gammas(2:num_species) + + ! Write fluxes vs Er + Write(str_num,*) num_ion_species + 2 ! Convert num to string + Write(iu_fvEr_out,'(f7.4,' // trim(adjustl(str_num)) // '(" ",e15.7))') & + roa_surf,Er_test/100._rknd,Gamma_e_vs_Er(ie),Gamma_i_vs_Er(ie,:) + + If ( output_QoT_vs_Er .EQV. .true. ) Then + QoT_e_vs_Er(ie) = QoTs(1) + QoT_i_vs_Er(ie,:) = QoTs(2:num_species) + Write(iu_QoTvEr_out,'(f7.4,' // trim(adjustl(str_num)) // '(" ",e15.7))') & + roa_surf,Er_test/100._rknd,QoT_e_vs_Er(ie),QoT_i_vs_Er(ie,:) + Endif + + ! Write flows vs Er + Write(str_num,*) (Smax+1)*num_species + 2 ! Convert num to string + Write(iu_flowvEr_out,'(f7.4,' // trim(adjustl(str_num)) // '(" ",e15.7))') & + roa_surf,Er_test/100._rknd,Flows + + Enddo !efield loop + RETURN + END SUBROUTINE penta_run_2_efield + + SUBROUTINE penta_run_3_ambipolar + USE vmec_var_pass + USE phys_const + USE coeff_var_pass + USE penta_functions_mod + USE PENTA_subroutines, ONLY: form_xvec, find_Er_roots + IMPLICIT NONE + ! Check for only one Er test value -- this is then used to evaluate the ambipolar fluxes QQ + !If ( num_Er_test == 1 ) Then + ! Er_roots = Er_test_vals + + ! Find the ambipolar root(s) from gamma_e = sum(Z*gamma_i) + Call find_Er_roots(gamma_e_vs_Er,gamma_i_vs_Er,Er_test_vals,Z_ion, & + num_Er_test,num_ion_species,Er_roots,num_roots) + + ! Allocate arrays according to number of ambipolar roots + Allocate(Flows_ambi((Smax+1)*num_species,num_roots)) ! Parallel flow moments + Allocate(Gammas_ambi(num_species,num_roots)) ! Rad. particle fluxes + Allocate(QoTs_ambi(num_species,num_roots)) ! Rad. energy fluxes + Allocate(Jprl_ambi(num_roots)) ! Parallel current density + Allocate(J_BS_ambi(num_roots)) ! BS current density + Allocate(sigma_par_ambi(num_roots)) ! Parallel conductivity + Allocate(sigma_par_Spitzer_ambi(num_roots)) ! Spitzer Parallel conductivity + Allocate(Jprl_parts(num_species,num_roots)) ! Par. curr. dens. per spec. + Allocate(upol(num_species,num_roots)) ! fsa contra pol flow + Allocate(utor(num_species,num_roots)) ! fsa contra tor flow + + ! Evaluate fluxes and flows at the ambipolar Er + Do iroot = 1_iknd, num_roots + + Er_test = Er_roots(iroot) + abs_Er = Dabs(Er_test) + + ! Form thermodynamic force vector (Xvec) + Call form_Xvec(Er_test,Z_ion,B_Eprl,num_ion_species,Xvec) + + ! Form alternate thermodynamic force vector (Avec) + Do ispec1 = 1,num_species + ind_X = (ispec1-1)*2 + 1 + ind_A = (ispec1-1)*3 + 1 + + Avec(ind_A) = -Xvec(ind_X) / (Temps(ispec1)*elem_charge) & + - 2.5_rknd*dTdrs(ispec1)/Temps(ispec1) + Avec(ind_A+1) = -Xvec(ind_X+1) / (Temps(ispec1)*elem_charge) + Avec(ind_A+2) = Xvec(num_species*2+1)*charges(ispec1) & + * B0/(Temps(ispec1)*elem_charge*Sqrt(Bsq)) + Enddo + + ! Select the appropriate algorithm and calculate the flows and fluxes + SELECT CASE (Method) + Case ('T', 'MBT') + ! Calculate array of parallel flow moments + ! Note: Flow methods are the same for T and MBT + Flows_ambi(:,iroot) = calc_flows_T(num_species,Smax,abs_Er,Temps,dens, & + vths,charges,masses,loglambda,B0,use_quanc8,Kmin,Kmax,numKsteps, & + log_interp,cmin,cmax,emin,emax,xt_c,xt_e,Dspl_D31,Dspl_logD33,num_c, & + num_e,kcord,keord,Avec,Bsq,lmat,J_BS) + ! Calculate array of radial particle fluxes + Gammas_ambi(:,iroot) = calc_fluxes_MBT(num_species,Smax,abs_Er,Temps, & + dens,vths,charges,masses,dTdrs,dndrs,loglambda,use_quanc8,Kmin,Kmax, & + numKsteps,log_interp,cmin,cmax,emin,emax,xt_c,xt_e,Dspl_logD11, & + Dspl_D31,Dspl_Dex,num_c,num_e,kcord,keord,Avec,lmat, & + Flows_ambi(:,iroot),U2,B0,flux_cap) + ! Calculate array of radial energy fluxes + QoTs_ambi(:,iroot) = calc_QoTs_MBT(num_species,Smax,abs_Er,Temps,dens, & + vths,charges,masses,dTdrs,dndrs,loglambda,use_quanc8,Kmin,Kmax, & + numKsteps,log_interp,cmin,cmax,emin,emax,xt_c,xt_e,Dspl_logD11, & + Dspl_D31,Dspl_Dex,num_c,num_e,kcord,keord,Avec,lmat, & + Flows_ambi(:,iroot),U2,B0,flux_cap) + + J_BS_ambi(iroot) = J_BS + Case ('SN') + ! Calculate array of parallel flow moments + + Flows_ambi(:,iroot) = calc_flows_SN(num_species,Smax,abs_Er,Temps,dens,& + vths,charges,masses,loglambda,B0,use_quanc8,Kmin,Kmax,numKsteps, & + log_interp,cmin,cmax,emin,emax,xt_c,xt_e,Dspl_Drat,Dspl_DUa,num_c, & + num_e,kcord,keord,Avec,lmat,sigma_par,sigma_par_Spitzer,J_BS) + ! Calculate array of radial particle fluxes + Gammas_ambi(:,iroot) = calc_fluxes_SN(num_species,Smax,abs_Er,Temps, & + dens,vths,charges,masses,loglambda,use_quanc8,Kmin,Kmax,numKsteps, & + log_interp,cmin,cmax,emin,emax,xt_c,xt_e,Dspl_Drat,Dspl_Drat2, & + Dspl_Dex,Dspl_logD11,Dspl_D31,num_c,num_e,kcord,keord,Avec,Bsq, & + lmat,Flows_ambi(:,iroot),U2,dTdrs,dndrs,flux_cap) + ! Calculate array of radial energy fluxes + QoTs_ambi(:,iroot) = calc_QoTs_SN(num_species,Smax,abs_Er,Temps,dens, & + vths,charges,masses,loglambda,use_quanc8,Kmin,Kmax,numKsteps, & + log_interp,cmin,cmax,emin,emax,xt_c,xt_e,Dspl_Drat,Dspl_Drat2, & + Dspl_Dex,Dspl_logD11,Dspl_D31,num_c,num_e,kcord,keord,Avec,Bsq, & + lmat,Flows_ambi(:,iroot),U2,dTdrs,dndrs,flux_cap) + + sigma_par_ambi(iroot) = sigma_par + sigma_par_Spitzer_ambi(iroot) = sigma_par_Spitzer + J_BS_ambi(iroot) = J_BS + Case ('DKES') + ! Calculate array of parallel flow moments + Flows_ambi(:,iroot) = calc_flows_DKES(num_species,Smax,abs_Er,Temps, & + dens,vths,charges,masses,loglambda,B0,use_quanc8,Kmin,Kmax,numKsteps,& + log_interp,cmin,cmax,emin,emax,xt_c,xt_e,Dspl_D31,Dspl_logD33,num_c, & + num_e,kcord,keord,Avec,J_BS) + ! Calculate array of radial particle fluxes + Gammas_ambi(:,iroot) = calc_fluxes_DKES(num_species,abs_Er,Temps,dens, & + vths,charges,masses,loglambda,use_quanc8,Kmin,Kmax,numKsteps, & + log_interp,cmin,cmax,emin,emax,xt_c,xt_e,Dspl_logD11,Dspl_D31,num_c, & + num_e,kcord,keord,Avec,B0) + ! Calculate array of radial energy fluxes + QoTs_ambi(:,iroot) = calc_QoTs_DKES(num_species,abs_Er,Temps,dens, & + vths,charges,masses,loglambda,use_quanc8,Kmin,Kmax,numKsteps, & + log_interp,cmin,cmax,emin,emax,xt_c,xt_e,Dspl_logD11,Dspl_D31,num_c, & + num_e,kcord,keord,Avec,B0) + + J_BS_ambi(iroot) = J_BS + Case Default + Write(*,'(3a)') ' Error: ''', Trim(Adjustl(Method)), & + ''' is not a valid Method' + Stop 'Error: Exiting, method select error in penta.f90 (4)' + ENDSELECT + + ! Calculate parallel current density + Jprl_parts(:,iroot) = dens*charges*Sqrt(Bsq) * & + Flows_ambi(1:(num_species-1)*(Smax+1)+1:Smax+1,iroot) + Jprl_ambi(iroot) = Sum(Jprl_parts(:,iroot)) + + ! Calculate flow components + upol(:,iroot) = chip*(4._rknd*pi*pi/vol_p)*( & + Flows_ambi(1:(num_species-1)*(Smax+1)+1:Smax+1,iroot) - & + bzeta*Xvec(1:(num_species-1)*2+1:2)/(charges*chip*Bsq)) + utor(:,iroot) = psip*(4._rknd*pi*pi/vol_p)*( & + Flows_ambi(1:(num_species-1)*(Smax+1)+1:Smax+1,iroot) + & + btheta*Xvec(1:(num_species-1)*2+1:2)/(charges*psip*Bsq)) + + END DO ! Ambipolar root loop + RETURN + END SUBROUTINE penta_run_3_ambipolar + + SUBROUTINE penta_run_4_cleanup + USE io_unit_spec + USE pprof_pass + USE vmec_var_pass + IMPLICIT NONE + ! First write output files + ! Loop over ambipolar Er for writing output files + Do iroot = 1_iknd, num_roots + + Er_test = Er_roots(iroot) + eaEr_o_kTe = arad*Er_test/Te + + ! Write fluxes to file "fluxes_vs_roa" + Write(str_num,*) 2*num_species + 2 + Write(iu_flux_out,'(f7.3,' // Trim(Adjustl(str_num)) // '(" ",e15.7))') & + roa_surf,Er_test/100._rknd,eaEr_o_kTe,Gammas_ambi(1,iroot), & + QoTs_ambi(1,iroot),Gammas_ambi(2:num_species,iroot), & + QoTs_ambi(2:num_species,iroot) + + ! Write flows to file "flows_vs_roa" + Write(str_num,*) (Smax+1)*num_species + 2 + Write(iu_flows_out,'(f7.3,' // trim(adjustl(str_num)) // '(" ",e15.7))') & + roa_surf,Er_test/100._rknd,eaEr_o_kTe,Flows_ambi(:,iroot) + + ! Write current densities to file "Jprl_vs_roa" + Write(str_num,*) num_species + 4 + Write(iu_Jprl_out,'(f7.3,' // trim(adjustl(str_num)) // '(" ",e15.7))') & + roa_surf,Er_test/100._rknd,eaEr_o_kTe,Jprl_parts(:,iroot),Jprl_ambi(iroot),J_BS_ambi(iroot) + + ! Write contravariant flows to file "ucontra_vs_roa" + Write(str_num,*) 2*num_species + 2 + Write(iu_contraflows_out,'(f7.3,' // trim(adjustl(str_num))//'(" ",e15.7))') & + roa_surf,Er_test/100._rknd,eaEr_o_kTe,upol(1,iroot),utor(1,iroot), & + upol(2:num_species,iroot),utor(2:num_species,iroot) + + ! Write sigmas to file "sigmas_vs_roa" + If( Method == 'SN') then + Write(str_num,*) 3 + Write(iu_sigmas_out,'(f7.3,' // trim(adjustl(str_num)) // '(" ",e15.7))') & + roa_surf,Er_test/100._rknd,sigma_par_ambi(iroot),sigma_par_Spitzer_ambi(iroot) + Endif + EndDo ! Ambipolar root loop + + ! Write plasma profile information to "plasma_profiles_check" + Write(str_num,*) 4*num_species + Write(iu_pprof_out,'(f7.3,' // trim(adjustl(str_num)) // '(" ",e15.7))') & + roa_surf,Te,ne,dnedr,dTedr,Ti,ni,dnidr,dTidr + + ! QQ write file with number of roots per surface! + + ! Write screen output + write(str_num,*) num_roots + write(*,'(f7.3,' // trim(adjustl(str_num)) // '(" ",e15.4))') & + roa_surf,er_roots(1:num_roots)/100._rknd + ! DEALLOCATE + CALL penta_deallocate_species + CALL penta_deallocate_dkescoeff + END SUBROUTINE penta_run_4_cleanup + + + +END MODULE PENTA_INTERFACE_MOD \ No newline at end of file From 9b8e4c5b6c8628be2536adfb68860c687d81c750 Mon Sep 17 00:00:00 2001 From: lazersos Date: Thu, 17 Oct 2024 10:29:43 +0200 Subject: [PATCH 19/47] PENTA: Change unit numbers from parameters to variables. --- PENTA/Sources/io_unit_spec.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PENTA/Sources/io_unit_spec.f90 b/PENTA/Sources/io_unit_spec.f90 index 375ee6115..061f5c107 100644 --- a/PENTA/Sources/io_unit_spec.f90 +++ b/PENTA/Sources/io_unit_spec.f90 @@ -21,7 +21,7 @@ Module io_unit_spec Implicit none - Integer, parameter :: & + Integer :: & iu_nl=21, & ! Ion parameter namelist file (input) iu_vmec=22, & ! VMEC data file (input) iu_pprof=23, & ! Plasma profile file (input) From cd12712f7e79d86cd5f8c22fa1d259b111eaf11d Mon Sep 17 00:00:00 2001 From: lazersos Date: Thu, 17 Oct 2024 13:46:49 +0200 Subject: [PATCH 20/47] PENTA: Subroutine interface checked. --- PENTA/Sources/penta_interface_mod.f90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/PENTA/Sources/penta_interface_mod.f90 b/PENTA/Sources/penta_interface_mod.f90 index ca80e2055..07b86cb6c 100644 --- a/PENTA/Sources/penta_interface_mod.f90 +++ b/PENTA/Sources/penta_interface_mod.f90 @@ -133,7 +133,7 @@ SUBROUTINE read_penta_run_params_namelist(filename,istat) END IF READ(iunit,NML=run_params,IOSTAT=istat) IF (istat /= 0) THEN - WRITE(6,'(A)') 'ERROR reading PENTA ion_params namelist from file: ',TRIM(filename) + WRITE(6,'(A)') 'ERROR reading PENTA run_params namelist from file: ',TRIM(filename) backspace(iunit) read(iunit,fmt='(A)') line write(6,'(A)') 'Invalid line in namelist: '//TRIM(line) @@ -280,7 +280,6 @@ SUBROUTINE penta_read_input_files USE coeff_var_pass, ONLY: D11_mat, cmul, num_c USE phys_const, ONLY: p_mass, elem_charge, e_mass USE read_input_file_mod - USE PENTA_subroutines, ONLY : define_friction_coeffs IMPLICIT NONE CALL read_vmec_file_2(js,run_ident) CALL read_pprof_file(pprof_char,num_ion_species,roa_surf,arad,kord_pprof) @@ -298,7 +297,7 @@ SUBROUTINE penta_read_input_files U2=1.5d0*D11_mat(num_c,1)/cmul(num_c); ENDIF ! Change Er test range to V/cm if necessary - IF ( input_is_Er) THEN + IF ( .not. input_is_Er) THEN Er_min = Er_min * Te / arad Er_max = Er_max * Te / arad ENDIF @@ -322,9 +321,6 @@ SUBROUTINE penta_read_input_files vths=(/vth_e,vth_i/) dTdrs=(/dTedr,dTidr/) dndrs=(/dnedr,dnidr/) - ! Define matrix of friction coefficients (lmat) - Call define_friction_coeffs(masses,charges,vths,Temps,dens,loglambda, & - num_species,Smax,lmat) RETURN END SUBROUTINE penta_read_input_files @@ -412,7 +408,7 @@ SUBROUTINE penta_open_output ! Set write status IF (i_append == 0) THEN fstatus = "unknown" - fpos = "asis" + fpos = "SEQUENTIAL" ELSEIF (i_append == 1) THEN fstatus = "old" fpos = "append" @@ -490,8 +486,11 @@ END SUBROUTINE penta_open_output SUBROUTINE penta_fit_rad_trans USE coeff_var_pass USE vmec_var_pass - USE PENTA_subroutines, ONLY : fit_coeffs + USE PENTA_subroutines, ONLY : fit_coeffs, define_friction_coeffs IMPLICIT NONE + ! Define matrix of friction coefficients (lmat) + Call define_friction_coeffs(masses,charges,vths,Temps,dens,loglambda, & + num_species,Smax,lmat) ! Fit radial transport coefficients specific to different methods SELECT CASE (Method) CASE ('T', 'MBT') @@ -537,6 +536,7 @@ END SUBROUTINE penta_fit_rad_trans SUBROUTINE penta_run_1_init IMPLICIT NONE INTEGER :: istat + CALL init_penta_input istat = 0 CALL read_penta_ion_params_namelist('ion_params',istat) istat = 0 From d94e7beba45768d33aec26f3094ec8b03f83af96 Mon Sep 17 00:00:00 2001 From: lazersos Date: Thu, 17 Oct 2024 14:23:17 +0200 Subject: [PATCH 21/47] PENTA: Added subroutines for outputting namelists. --- PENTA/Sources/penta_interface_mod.f90 | 91 +++++++++++++++++++++++++++ 1 file changed, 91 insertions(+) diff --git a/PENTA/Sources/penta_interface_mod.f90 b/PENTA/Sources/penta_interface_mod.f90 index 07b86cb6c..24480b1e1 100644 --- a/PENTA/Sources/penta_interface_mod.f90 +++ b/PENTA/Sources/penta_interface_mod.f90 @@ -114,6 +114,43 @@ SUBROUTINE read_penta_ion_params_namelist(filename,istat) RETURN END SUBROUTINE read_penta_ion_params_namelist + SUBROUTINE write_ion_params_nml(iunit) + IMPLICIT NONE + INTEGER(iknd), INTENT(inout) :: iunit + CHARACTER(LEN=*), PARAMETER :: outint = "(2X,A,1X,'=',1X,I0)" + INTEGER(iknd) :: k + WRITE(iunit,'(A)') '&ION_PARAMS' + WRITE(iunit,'(A)') '----------------------------------------------------------------' + WRITE(iunit,outint) 'NUM_ION_SPECIES',num_ion_species + WRITE(iunit,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'Z_ION_INIT',(Z_ion_init(k), k=1,num_ion_species) + WRITE(iunit,"(2X,A,1X,'=',4(1X,ES22.12E3))") 'MIOMP_INIT',(miomp_init(k), k=1,num_ion_species) + WRITE(iunit,'(A)') '/\n' + END SUBROUTINE write_ion_params_nml + + SUBROUTINE write_ion_params_namelist_byfile(filename) + USE safe_open_mod + IMPLICIT NONE + CHARACTER(LEN=*), INTENT(in) :: filename + INTEGER(iknd) :: iunit, istat + LOGICAL :: lexists + + iunit = 100 + istat = 0 + INQUIRE(FILE=TRIM(filename),exist=lexists) + IF (lexists) THEN + CALL safe_open(iunit,istat,TRIM(filename),'old','formatted', access_in='append') + !OPEN(unit=iunit, file=TRIM(filename),iostat=istat, status="old", position="append") + ELSE + CALL safe_open(iunit,istat,TRIM(filename),'replace','formatted') + !OPEN(unit=iunit, file=TRIM(filename),iostat=istat, status="new") + END IF + IF (istat .ne. 0) RETURN + CALL write_ion_params_nml(iunit) + CLOSE(iunit) + + RETURN + END SUBROUTINE write_ion_params_namelist_byfile + SUBROUTINE read_penta_run_params_namelist(filename,istat) USE safe_open_mod IMPLICIT NONE @@ -144,6 +181,60 @@ SUBROUTINE read_penta_run_params_namelist(filename,istat) RETURN END SUBROUTINE read_penta_run_params_namelist + SUBROUTINE write_run_params_nml(iunit) + IMPLICIT NONE + INTEGER(iknd), INTENT(inout) :: iunit + CHARACTER(LEN=*), PARAMETER :: outboo = "(2X,A,1X,'=',1X,L1)" + CHARACTER(LEN=*), PARAMETER :: outint = "(2X,A,1X,'=',1X,I0)" + CHARACTER(LEN=*), PARAMETER :: outdbl = "(2X,A,1X,'=',1X,ES22.14)" + CHARACTER(LEN=*), PARAMETER :: outstr = "(2X,A,1X,'=',1X,'''',A,'''')" + INTEGER(iknd) :: k + WRITE(iunit,'(A)') '&RUN_PARAMS' + WRITE(iunit,outboo) 'INPUT_IS_ER',input_is_er + WRITE(iunit,outboo) 'LOG_INTERP',log_interp + WRITE(iunit,outboo) 'USE_QUANC8',use_quanc8 + WRITE(iunit,outboo) 'READ_U2_FILE',read_U2_file + WRITE(iunit,outboo) 'ADD_SPITZER_TO_D33',Add_Spitzer_to_D33 + WRITE(iunit,outboo) 'FLUX_CAP',flux_cap + WRITE(iunit,outboo) 'OUTPUT_QOT_VS_ER',output_QoT_vs_Er + WRITE(iunit,outboo) 'USE_BEAM',use_beam + WRITE(iunit,outint) 'NUM_ER_TEST',num_Er_test + WRITE(iunit,outint) 'NUMKSTEPS',numksteps + WRITE(iunit,outint) 'KORD_PPROF',kord_pprof + WRITE(iunit,outint) 'KEORD',keord + WRITE(iunit,outint) 'KCORD',kcord + WRITE(iunit,outdbl) 'KMIN',kmin + WRITE(iunit,outdbl) 'KMAX',kmax + WRITE(iunit,outdbl) 'EPSABS',epsabs + WRITE(iunit,outdbl) 'EPSREL',epsrel + WRITE(iunit,outstr) 'METHOD',method + WRITE(iunit,'(A)') '/\n' + END SUBROUTINE write_run_params_nml + + SUBROUTINE write_run_params_namelist_byfile(filename) + USE safe_open_mod + IMPLICIT NONE + CHARACTER(LEN=*), INTENT(in) :: filename + INTEGER(iknd) :: iunit, istat + LOGICAL :: lexists + + iunit = 100 + istat = 0 + INQUIRE(FILE=TRIM(filename),exist=lexists) + IF (lexists) THEN + CALL safe_open(iunit,istat,TRIM(filename),'old','formatted', access_in='append') + !OPEN(unit=iunit, file=TRIM(filename),iostat=istat, status="old", position="append") + ELSE + CALL safe_open(iunit,istat,TRIM(filename),'replace','formatted') + !OPEN(unit=iunit, file=TRIM(filename),iostat=istat, status="new") + END IF + IF (istat .ne. 0) RETURN + CALL write_run_params_nml(iunit) + CLOSE(iunit) + + RETURN + END SUBROUTINE write_run_params_namelist_byfile + SUBROUTINE penta_init_commandline IMPLICIT NONE CALL GETCARG(1, arg1, numargs) From 61a523d36ee8b63a936e2e575901b4b5372cf4bb Mon Sep 17 00:00:00 2001 From: lazersos Date: Thu, 17 Oct 2024 14:35:15 +0200 Subject: [PATCH 22/47] PENTA: Removed extra files. --- PENTA/Sources/penta_input_mod.f | 103 ------------------------ PENTA/Sources/penta_main.f | 135 -------------------------------- 2 files changed, 238 deletions(-) delete mode 100644 PENTA/Sources/penta_input_mod.f delete mode 100644 PENTA/Sources/penta_main.f diff --git a/PENTA/Sources/penta_input_mod.f b/PENTA/Sources/penta_input_mod.f deleted file mode 100644 index e4640084d..000000000 --- a/PENTA/Sources/penta_input_mod.f +++ /dev/null @@ -1,103 +0,0 @@ - module penta_input_mod - use penta_kind_mod - implicit none - INTEGER(iknd), parameter :: num_ion_max = 20_iknd !Maximum number of ion species to be read from namelist file - INTEGER(iknd) :: num_ion_species ! Number of ion species - REAL(rknd), DIMENSION(num_ion_max) :: Z_ion_init,miomp_init - - namelist / ion_params / num_ion_species,Z_ion_init,miomp_init - contains - - SUBROUTINE init_ion_params_nml - IMPLICIT NONE - num_ion_species = 1 - Z_ion_init = 1.0_rknd - miomp_init = 1.0_rknd - END SUBROUTINE init_ion_params_nml - - SUBROUTINE read_ion_params_nml(filename,istat) - USE safe_open_mod - IMPLICIT NONE - ! Subroutine parameters - CHARACTER(256), INTENT(INOUT) :: filename - INTEGER(iknd), INTENT(INOUT) :: istat - ! Local Variables - LOGICAL :: lexist - INTEGER(iknd) :: iunit - CHARACTER(256) :: filename_local - - ! Handle filename - filename_local = TRIM(filename) - - !Check which file is available - INQUIRE(FILE=TRIM(filename_local),EXIST=lexist) - IF (.not. lexist) THEN - WRITE(6,*) ' -Could not open input file trying ion_params' - filename_local = 'ion_params' - INQUIRE(FILE=TRIM(filename_local),EXIST=lexist) - IF (.not. lexist) THEN - WRITE(6,*) ' -Could not open ion_params file' - istat = -1 - RETURN - ENDIF - ENDIF - - ! Open and read the file - CALL safe_open(iunit,istat,TRIM(filename_local),'old', - 1 'formatted') - IF (istat /= 0) THEN - WRITE(6,*) 'ERROR: Could no open: ',TRIM(filename_local) - RETURN - ENDIF - READ(iunit,NML=ion_params,IOSTAT=istat) - IF (istat /= 0) THEN - WRITE(6,*) 'ERROR: ion_params in ',TRIM(filename_local) - !backspace(iu_nl) - !read(iu_nl,fmt='(A)') line - !write(6,'(A)') 'Invalid line in namelist: '//TRIM(line) - RETURN - ENDIF - CLOSE(iunit) - RETURN - END SUBROUTINE read_ion_params_nml - - SUBROUTINE write_ion_params_nml(iunit) - IMPLICIT NONE - INTEGER(iknd), INTENT(inout) :: iunit - CHARACTER(LEN=*), PARAMETER :: outint = "(2X,A,1X,'=',1X,I0)" - INTEGER(iknd) :: k - WRITE(iunit,'(A)') '&ION_PARAMS' - WRITE(iunit,'(A)') - 1'----------------------------------------------------------------' - WRITE(iunit,outint) 'NUM_ION_SPECIES',num_ion_species - WRITE(iunit,"(2X,A,1X,'=',4(1X,ES22.12E3))") - 1'Z_ION_INIT',(Z_ion_init(k), k=1,num_ion_species) - WRITE(iunit,"(2X,A,1X,'=',4(1X,ES22.12E3))") - 1'MIOMP_INIT',(miomp_init(k), k=1,num_ion_species) - WRITE(iunit,'(A)') '/\n' - END SUBROUTINE write_ion_params_nml - - SUBROUTINE write_ion_params_namelist_byfile(filename) - USE safe_open_mod - IMPLICIT NONE - CHARACTER(LEN=*), INTENT(in) :: filename - INTEGER(iknd) :: iunit, istat - LOGICAL :: lexists - - iunit = 100 - istat = 0 - INQUIRE(FILE=TRIM(filename),exist=lexists) - IF (lexists) THEN - OPEN(unit=iunit, file=TRIM(filename), - 1 iostat=istat, status="old", position="append") - ELSE - OPEN(unit=iunit, file=TRIM(filename), - 1 iostat=istat, status="new") - END IF - IF (istat .ne. 0) RETURN - CALL write_ion_params_nml(iunit) - CLOSE(iunit) - - RETURN - END SUBROUTINE write_ion_params_namelist_byfile - end module penta_input_mod \ No newline at end of file diff --git a/PENTA/Sources/penta_main.f b/PENTA/Sources/penta_main.f deleted file mode 100644 index e49aa383a..000000000 --- a/PENTA/Sources/penta_main.f +++ /dev/null @@ -1,135 +0,0 @@ -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c -c -PENTA with Impurities- -c -ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c To run,type: PENTA_imp [command line arguments] -c -c where the command line arguments are: -c -c 1 unique identifier text on DKES database files -c i.e., for (m,l,n)star_lijs_***, this would be -c the text that replaces *** for the specific surface (ex: hsx_s10) -c 2 Er,min -c 3 Er,max - args 2 and 3 specify the interval in Er over which -c the search is done for the ambipolar electric field root. -c Er,min and Er,max are dimensionless normalized electric -c field variables, i.e.: eEr,min/max/kT_e -c - Note if the logical variable input_is_Er below is true then -c the search interval is set by args 2 and 3 as Er (V/cm) -c 4 flux surface number -c 5 parameter that determines whether output data should be -c written into new files (=0) or appended to -c existing files (=1). Typically,when a script is run for a -c sequence of flux surfaces, this is set to 0 for the first -c surface and to 1 for subsequent surfaces -c 6 extension of the profile_data_*** file (i.e., the *** text) -c that comes from extracting the appropriate profiles from -c the VMEC run using the profile_extractor.f code -c 7 unique identifier on plasma profiles file, i.e., the *** -c text in the filename plasma_profiles***.dat. This allows -c multiple profiles to be run without having to rename files -c 8 in T*V/m -c -c -c Required input files: -c -c [l,m,n]star_lijs_***_s# - where *** is arg1 above and # arg 4. -c These files contain the values of efield and cmul and the -c corresponding normalized monoenergetic coefficients. Note -c that L* should be L*/e**2 where e is the elementary charge. -c -c ion_params - A file containing the namelist "ion_params" which -c defines the variables num_ion_species, Z_ion_init and -c miomp_init which are the number of ion species (integer), -c the corresponding ion charge numbers (real, array) and -c the ion to proton mass ratio (real, array) for each -c non-electron species. -c -c profile_data_*** - where *** is arg1 above. Contains geometry -c variables, most likely from a VMEC run. -c -c plasma_profilesXXX.dat - A file containing the plasma profile -c data, where XXX is arg8 above. This file has the format: -c row 1: number of radial points for which data is specified. -c all other rows: r/a, ne, Te, Ti1, ni1, Ti2, ni2 ... -c where the densities are in units of 10**12/cc and the -c temperatures in eV and i1, i2... are the ion species. -c -c Utilde2_profile - Contains the quantity , where U is the -c PS flow function as defined in Sugama and Nishimura. The -c first row is the number of points, then r/a points and the -c corresponding value. -c -c Output files: -c -c "fluxes_vs_roa", "fluxes_vs_Er", "flows_vs_roa", "plasma_profiles_check" -c -c -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c NOTES -c - Be sure to check parameters below -c - based on Sugama, Nishimura PoP 9 4637 (2002) and -c Sugama, Nishimura PoP 15, 042502 (2008) -c - Some code used from orignal PENTA code, written by Don Spong and -c modified by Jeremy Lore -c -c 7/2009 JL -ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c PROGRAMMING NOTES -c -speed optimization -c -recursive ambipolar solver -c -fix numargs -ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -c 2024 Notes -c - Created PENTA_IMP_MAIN program and changed PENTA_IMP to a -c callable subroutine. -c - Can read wout files instead of profile_data_**** -c - Routine to set profiles instead of reading plasma_profiles_**** -c - Routine to set L/M/NIJ coefficients instead of reading them. -c - Routine to set Utilde2 value instead of reading from file. -c -cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - PROGRAM PENTA_MAIN - use penta_kind_mod - USE read_wout_mod, ONLY: read_wout_deallocate - IMPLICIT NONE - integer(iknd) :: numargs, js, i_append - real(rknd) :: Er_min, Er_max, B_Eprl - character(100) :: arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8 - character(100) :: coeff_ext, run_ident - character(10) :: pprof_char - numargs = 8 - call getarg(1, arg1) - call getarg(2, arg2) - call getarg(3, arg3) - call getarg(4, arg4) - call getarg(5, arg5) - call getarg(6, arg6) - call getarg(7, arg7) - if (numargs .eq. 8) call getarg(8, arg8) - - !store command line args - coeff_ext = trim(adjustl(arg1)) - read(arg2,*) Er_min - read(arg3,*) Er_max - read(arg4,*) js - read(arg5,*) i_append - run_ident = trim(adjustl(arg6)) - pprof_char = trim(adjustl(arg7)) - if (numargs .eq. 8) read(arg8,*) B_Eprl - - ! Call subroutine - CALL PENTA_IMP(coeff_ext, Er_min, Er_max, js, i_append, - 1 run_ident, pprof_char, B_Eprl) - - ! DEALLOCATE VMEC DATA - CALL read_wout_deallocate - - END PROGRAM PENTA_MAIN - From 1e750f7548608105e611c24b8d448cd550bff7f7 Mon Sep 17 00:00:00 2001 From: lazersos Date: Thu, 17 Oct 2024 15:02:56 +0200 Subject: [PATCH 23/47] THIFT: Added prototype of call to PENTA. --- THRIFT/Sources/thrift_penta.f90 | 77 +++++++++++++++++++++++++++++++++ 1 file changed, 77 insertions(+) create mode 100644 THRIFT/Sources/thrift_penta.f90 diff --git a/THRIFT/Sources/thrift_penta.f90 b/THRIFT/Sources/thrift_penta.f90 new file mode 100644 index 000000000..15786be24 --- /dev/null +++ b/THRIFT/Sources/thrift_penta.f90 @@ -0,0 +1,77 @@ +!----------------------------------------------------------------------- +! Subroutine: thrift_penta +! Authors: S. Lazerson (samuel.lazerson@gauss-fusion.com) +! Date: 08/22/2024 +! Description: This subroutine calculates the PENTA data. +!----------------------------------------------------------------------- + SUBROUTINE thrift_penta(lscreen,iflag) +!----------------------------------------------------------------------- +! Libraries +!----------------------------------------------------------------------- + USE thrift_runtime + USE thrift_input_mod + USE thrift_vars, nrho_thrift => nrho + USE thrift_profiles_mod + USE thrift_equil + USE thrift_funcs + USE penta_interface_mod + USE mpi_params + USE mpi_inc +!----------------------------------------------------------------------- +! Subroutine Parameters +! lscreen Screen output +! iflag Error flag +!----------------------------------------------------------------------- + IMPLICIT NONE + LOGICAL, INTENT(in) :: lscreen + INTEGER, INTENT(inout) :: iflag +!----------------------------------------------------------------------- +! Local Variables +!----------------------------------------------------------------------- + INTEGER :: ns_dkes, k + REAL(rprec) :: s +!----------------------------------------------------------------------- +! BEGIN SUBROUTINE +!----------------------------------------------------------------------- + IMPLICIT NONE + IF (iflag < 0) RETURN + IF (lscreen) WRITE(6,'(a)') ' -------------------- NEOCLASSICAL FLUX USING PENTA -------------------' + IF (lvmec) THEN + ierr_mpi = 0 + ! PENTA is parallelized over radial surfaces in this routine. + ns_dkes = 0 + DO k = 1, DKES_NS_MAX + IF ((DKES_K(k) > 0)) ns_dkes = ns_dkes+1 + END DO + ! Break up work + CALL MPI_CALC_MYRANGE(MPI_COMM_MYWORLD,1,ns_dkes,mystart,myend) + DO k = mystart,myend + ! This part mimics penta_run_1 but without reading the namelists + Er_min = -250.0 + Er_max = 250.0 + js = DKES_K(k) + i_append = 0 + B_Eprl = 0.0 + Smax = 1 + coeff_ext = '' + run_ident = '' + pprof_char = '' + CALL PENTA_ALLOCATE_SPECIES + CALL PENTA_READ_INPUT_FILES + CALL PENTA_SCREEN_INFO + CALL PENTA_ALLOCATE_DKESCOEFF + CALL PENTA_FIT_DXX_COEF + CALL PENTA_OPEN_OUTPUT + CALL PENTA_FIT_RAD_TRANS + ! Now the basic steps + CALL PENTA_RUN_2_EFIELD + CALL PENTA_RUN_3_AMBIPOLAR + CALL PENTA_RUN_4_CLEANUP + END DO + ENDIF + IF (lscreen) WRITE(6,'(a)') ' ------------------- NEOCLASSICAL FLUX CALCULATION DONE ---------------------' + RETURN +!----------------------------------------------------------------------- +! END SUBROUTINE +!----------------------------------------------------------------------- + END SUBROUTINE thrift_penta \ No newline at end of file From 515faedbabc40b24faa3e0ead2db74cb06e820dc Mon Sep 17 00:00:00 2001 From: lazersos Date: Thu, 17 Oct 2024 15:41:40 +0200 Subject: [PATCH 24/47] PENTA: Added unique output file names possible. --- PENTA/Sources/penta_interface_mod.f90 | 28 +++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/PENTA/Sources/penta_interface_mod.f90 b/PENTA/Sources/penta_interface_mod.f90 index 24480b1e1..f9f97a1df 100644 --- a/PENTA/Sources/penta_interface_mod.f90 +++ b/PENTA/Sources/penta_interface_mod.f90 @@ -490,11 +490,13 @@ SUBROUTINE penta_screen_info RETURN END SUBROUTINE penta_screen_info - SUBROUTINE penta_open_output + SUBROUTINE penta_open_output(file_ext) USE safe_open_mod USE io_unit_spec IMPLICIT NONE + CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: file_ext INTEGER :: istat + CHARACTER(LEN=256) :: local_ext istat = 0 ! Set write status IF (i_append == 0) THEN @@ -508,35 +510,41 @@ SUBROUTINE penta_open_output CALL FLUSH(6) STOP 'Error: Exiting, i_append error in penta.f90' END IF + ! Create local extension + IF (PRESENT(file_ext)) THEN + local_ext = "."//TRIM(file_ext) + ELSE + local_ext = '' + END IF ! Open files !Open(unit=iu_flux_out, file="fluxes_vs_roa", position=Trim(Adjustl(fpos)),status=Trim(Adjustl(fstatus))) - CALL safe_open(iu_flux_out, istat, "fluxes_vs_roa", & + CALL safe_open(iu_flux_out, istat, "fluxes_vs_roa"//TRIM(local_ext), & Trim(Adjustl(fstatus)), 'formatted',& access_in=Trim(Adjustl(fpos))) - CALL safe_open(iu_pprof_out, istat, "plasma_profiles_check", & + CALL safe_open(iu_pprof_out, istat, "plasma_profiles_check"//TRIM(local_ext), & Trim(Adjustl(fstatus)), 'formatted',& access_in=Trim(Adjustl(fpos))) - CALL safe_open(iu_fvEr_out, istat, "fluxes_vs_Er", & + CALL safe_open(iu_fvEr_out, istat, "fluxes_vs_Er"//TRIM(local_ext), & Trim(Adjustl(fstatus)), 'formatted',& access_in=Trim(Adjustl(fpos))) - CALL safe_open(iu_flows_out, istat, "flows_vs_roa", & + CALL safe_open(iu_flows_out, istat, "flows_vs_roa"//TRIM(local_ext), & Trim(Adjustl(fstatus)), 'formatted',& access_in=Trim(Adjustl(fpos))) - CALL safe_open(iu_flowvEr_out, istat, "flows_vs_Er", & + CALL safe_open(iu_flowvEr_out, istat, "flows_vs_Er"//TRIM(local_ext), & Trim(Adjustl(fstatus)), 'formatted',& access_in=Trim(Adjustl(fpos))) - CALL safe_open(iu_Jprl_out, istat, "Jprl_vs_roa", & + CALL safe_open(iu_Jprl_out, istat, "Jprl_vs_roa"//TRIM(local_ext), & Trim(Adjustl(fstatus)), 'formatted',& access_in=Trim(Adjustl(fpos))) - CALL safe_open(iu_contraflows_out, istat, "ucontra_vs_roa", & + CALL safe_open(iu_contraflows_out, istat, "ucontra_vs_roa"//TRIM(local_ext), & Trim(Adjustl(fstatus)), 'formatted',& access_in=Trim(Adjustl(fpos))) IF (method == 'SN') & - CALL safe_open(iu_sigmas_out, istat, "sigmas_vs_roa", & + CALL safe_open(iu_sigmas_out, istat, "sigmas_vs_roa"//TRIM(local_ext), & Trim(Adjustl(fstatus)), 'formatted',& access_in=Trim(Adjustl(fpos))) IF (output_QoT_vs_Er) THEN - CALL safe_open(iu_QoTvEr_out, istat, "QoTs_vs_Er", & + CALL safe_open(iu_QoTvEr_out, istat, "QoTs_vs_Er"//TRIM(local_ext), & Trim(Adjustl(fstatus)), 'formatted',& access_in=Trim(Adjustl(fpos))) Write(iu_QoTvEr_out,'("*",/,"r/a Er[V/cm] Q_e/T_e [m**-2s**-1] ",& From 0f10b84e6b5f04ddc812c15fd1d6442427477a78 Mon Sep 17 00:00:00 2001 From: lazersos Date: Thu, 17 Oct 2024 15:42:05 +0200 Subject: [PATCH 25/47] THRIFT: PENTA interface uses unique output file names. --- THRIFT/Sources/thrift_penta.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/THRIFT/Sources/thrift_penta.f90 b/THRIFT/Sources/thrift_penta.f90 index 15786be24..e63620b83 100644 --- a/THRIFT/Sources/thrift_penta.f90 +++ b/THRIFT/Sources/thrift_penta.f90 @@ -50,7 +50,7 @@ SUBROUTINE thrift_penta(lscreen,iflag) Er_min = -250.0 Er_max = 250.0 js = DKES_K(k) - i_append = 0 + i_append = 1 B_Eprl = 0.0 Smax = 1 coeff_ext = '' @@ -61,7 +61,7 @@ SUBROUTINE thrift_penta(lscreen,iflag) CALL PENTA_SCREEN_INFO CALL PENTA_ALLOCATE_DKESCOEFF CALL PENTA_FIT_DXX_COEF - CALL PENTA_OPEN_OUTPUT + CALL PENTA_OPEN_OUTPUT(proc_string) CALL PENTA_FIT_RAD_TRANS ! Now the basic steps CALL PENTA_RUN_2_EFIELD From 94330d2b2c598f14c57e0417fb62fa2a60d28096 Mon Sep 17 00:00:00 2001 From: lazersos Date: Fri, 18 Oct 2024 12:29:31 +0200 Subject: [PATCH 26/47] PENTA: Added wrapper routines for ion namelist, command line, and vmec data --- PENTA/Sources/penta_interface_mod.f90 | 81 +++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) diff --git a/PENTA/Sources/penta_interface_mod.f90 b/PENTA/Sources/penta_interface_mod.f90 index f9f97a1df..b286ed274 100644 --- a/PENTA/Sources/penta_interface_mod.f90 +++ b/PENTA/Sources/penta_interface_mod.f90 @@ -82,6 +82,17 @@ SUBROUTINE init_penta_input RETURN END SUBROUTINE init_penta_input + SUBROUTINE penta_set_ion_params(num_ion_in, Z_ion_in, miomp_in) + IMPLICIT NONE + INTEGER, INTENT(in) :: num_ion_in + REAL(rknd), DIMENSION(NUM_ION_MAX), INTENT(in) :: Z_ion_in + REAL(rknd), DIMENSION(NUM_ION_MAX), INTENT(in) :: miomp_in + num_ion_species = num_ion_in + Z_ion_init = Z_ion_in + miomp_init = miomp_in + RETURN + END SUBROUTINE penta_set_ion_params + SUBROUTINE read_penta_ion_params_namelist(filename,istat) USE safe_open_mod IMPLICIT NONE @@ -264,6 +275,29 @@ SUBROUTINE penta_init_commandline RETURN END SUBROUTINE penta_init_commandline + SUBROUTINE penta_set_commandline(Er_min_in,Er_max_in,js_in,i_append_in,B_Eprl_in,Smax_in,ext_in,run_in,pprof_in) + IMPLICIT NONE + REAL(rknd), INTENT(IN) :: Er_min_in + REAL(rknd), INTENT(IN) :: Er_max_in + INTEGER(iknd), INTENT(IN) :: js_in + INTEGER(iknd), INTENT(IN) :: i_append_in + REAL(rknd), INTENT(IN) :: B_Eprl_in + REAL(rknd), INTENT(IN) :: Smax_in + CHARACTER(LEN=*), INTENT(IN) :: ext_in + CHARACTER(LEN=*), INTENT(IN) :: run_in + CHARACTER(LEN=*), INTENT(IN) :: pprof_in + Er_min = Er_min_in + Er_max = Er_max_in + js = js_in + i_append = i_append_in + B_Eprl = B_Eprl_in + Smax = Smax_in + coeff_ext = TRIM(ADJUSTL(ext_in)) + run_ident = TRIM(ADJUSTL(run_in)) + pprof_char = TRIM(ADJUSTL(pprof_in)) + RETURN + END SUBROUTINE penta_set_commandline + SUBROUTINE penta_allocate_species USE pprof_pass, ONLY: ni, Ti, dnidr, dTidr IMPLICIT NONE @@ -415,6 +449,53 @@ SUBROUTINE penta_read_input_files RETURN END SUBROUTINE penta_read_input_files + SUBROUTINE penta_set_eq_data(rho,A_in,R_in,V_in,chip_in,phip_in,iota_in,bth_in,bz_in,Bsq_in,B0_in) + USE vmec_var_pass + IMPLICIT NONE + REAL(rknd), INTENT(IN) :: rho + REAL(rknd), INTENT(IN) :: A_in + REAL(rknd), INTENT(IN) :: R_in + REAL(rknd), INTENT(IN) :: V_in + REAL(rknd), INTENT(IN) :: chip_in + REAL(rknd), INTENT(IN) :: phip_in + REAL(rknd), INTENT(IN) :: iota_in + REAL(rknd), INTENT(IN) :: bth_in + REAL(rknd), INTENT(IN) :: bze_in + REAL(rknd), INTENT(IN) :: Bsq_in + REAL(rknd), INTENT(IN) :: B0_in + REAL(rknd) :: TWOPI + TWOPI = 8*ATAN(1.0_rknd) + roa_surf = rho + arad = A_in + Rmajor = R_in + r_surf = roa_surf*arad + chip = chip_in + psip = phip_in + iota = iota_in + btheta = bth_in + bzeta = bze_in + Bsq = Bsq_in + vol_p = V_in + B0 = B0_in + ! Note that vp from VMEC comes normalized by 4pi^2 + ! Therefore we need to denormalize it + vol_p = TWOPI*TWOPI*vol_p + ! vp_vmec is ~dV/ds, but what we want is dVdr=dVds*dsdr + ! Since PENTA uses r/a=sqrt(s), then dVdr=dVds*2*r/a^2 + vol_p = vol_p * 2.0_rknd*r_surf/arad**2 + !Same for psip and chip: need to convert d/ds to d/dr + psip = psip * 2.0_rknd*r_surf/arad**2 + chip = chip * 2.0_rknd*r_surf/arad**2 + !psip and chip are used to compute flows in theta and zeta direction (see penta.f90) + !It is assumed that psip and chip are normalized by 2pi + !So need to do it here: + psip = psip / TWOPI + chip = chip / TWOPI + ! For the SN formulation this is the think to do, but not for the the other formulations... + b0 = dsqrt(bsq) + RETURN + END SUBROUTINE penta_set_eq_data + SUBROUTINE penta_fit_DXX_coef USE coeff_var_pass USE PENTA_subroutines, ONLY : fit_coeffs From 6608bc53b1ef32557aef35f71519e632d674f902 Mon Sep 17 00:00:00 2001 From: lazersos Date: Fri, 18 Oct 2024 14:16:25 +0200 Subject: [PATCH 27/47] PENTA: Added the missing interfaces to avoid using files. --- PENTA/Sources/penta_interface_mod.f90 | 125 ++++++++++++++++++++++---- 1 file changed, 108 insertions(+), 17 deletions(-) diff --git a/PENTA/Sources/penta_interface_mod.f90 b/PENTA/Sources/penta_interface_mod.f90 index b286ed274..702d40be5 100644 --- a/PENTA/Sources/penta_interface_mod.f90 +++ b/PENTA/Sources/penta_interface_mod.f90 @@ -399,27 +399,36 @@ SUBROUTINE penta_deallocate_dkescoeff RETURN END SUBROUTINE penta_deallocate_dkescoeff - SUBROUTINE penta_read_input_files + SUBROUTINE penta_read_input_files(lvmec,lpprof,ldkes,lbeam,lU2) USE vmec_var_pass, ONLY: roa_surf, arad, Bsq USE pprof_pass USE coeff_var_pass, ONLY: D11_mat, cmul, num_c USE phys_const, ONLY: p_mass, elem_charge, e_mass USE read_input_file_mod IMPLICIT NONE - CALL read_vmec_file_2(js,run_ident) - CALL read_pprof_file(pprof_char,num_ion_species,roa_surf,arad,kord_pprof) - CALL read_dkes_star_files(coeff_ext,Add_Spitzer_to_D33,Bsq) - IF (use_beam) THEN - CALL read_beam_file(roa_surf,kord_pprof) - ELSE - beam_force = 0._rknd + LOGICAL, INTENT(IN) :: lvmec + LOGICAL, INTENT(IN) :: lpprof + LOGICAL, INTENT(IN) :: ldkes + LOGICAL, INTENT(IN) :: lbeam + LOGICAL, INTENT(IN) :: lU2 + IF (lvmec) CALL read_vmec_file_2(js,run_ident) + IF (lpprof) CALL read_pprof_file(pprof_char,num_ion_species,roa_surf,arad,kord_pprof) + IF (ldkes) CALL read_dkes_star_files(coeff_ext,Add_Spitzer_to_D33,Bsq) + IF (lbeam) THEN + IF (use_beam) THEN + CALL read_beam_file(roa_surf,kord_pprof) + ELSE + beam_force = 0._rknd + ENDIF ENDIF ! Optionally read file containing info. Else this is ! calculated from the D11* coefficient at high nu/v and Er=0. - IF ( read_U2_file ) THEN - CALL read_Utilde2_file(roa_surf,U2,kord_pprof) - ELSE - U2=1.5d0*D11_mat(num_c,1)/cmul(num_c); + IF (lU2) THEN + IF ( read_U2_file) THEN + CALL read_Utilde2_file(roa_surf,U2,kord_pprof) + ELSE + U2=1.5d0*D11_mat(num_c,1)/cmul(num_c); + ENDIF ENDIF ! Change Er test range to V/cm if necessary IF ( .not. input_is_Er) THEN @@ -449,7 +458,7 @@ SUBROUTINE penta_read_input_files RETURN END SUBROUTINE penta_read_input_files - SUBROUTINE penta_set_eq_data(rho,A_in,R_in,V_in,chip_in,phip_in,iota_in,bth_in,bz_in,Bsq_in,B0_in) + SUBROUTINE penta_set_eq_data(rho,A_in,R_in,V_in,chip_in,phip_in,iota_in,bth_in,bze_in,Bsq_in) USE vmec_var_pass IMPLICIT NONE REAL(rknd), INTENT(IN) :: rho @@ -462,7 +471,6 @@ SUBROUTINE penta_set_eq_data(rho,A_in,R_in,V_in,chip_in,phip_in,iota_in,bth_in,b REAL(rknd), INTENT(IN) :: bth_in REAL(rknd), INTENT(IN) :: bze_in REAL(rknd), INTENT(IN) :: Bsq_in - REAL(rknd), INTENT(IN) :: B0_in REAL(rknd) :: TWOPI TWOPI = 8*ATAN(1.0_rknd) roa_surf = rho @@ -476,7 +484,6 @@ SUBROUTINE penta_set_eq_data(rho,A_in,R_in,V_in,chip_in,phip_in,iota_in,bth_in,b bzeta = bze_in Bsq = Bsq_in vol_p = V_in - B0 = B0_in ! Note that vp from VMEC comes normalized by 4pi^2 ! Therefore we need to denormalize it vol_p = TWOPI*TWOPI*vol_p @@ -492,10 +499,94 @@ SUBROUTINE penta_set_eq_data(rho,A_in,R_in,V_in,chip_in,phip_in,iota_in,bth_in,b psip = psip / TWOPI chip = chip / TWOPI ! For the SN formulation this is the think to do, but not for the the other formulations... - b0 = dsqrt(bsq) + b0 = dsqrt(Bsq) RETURN END SUBROUTINE penta_set_eq_data + SUBROUTINE penta_set_pprof(ne_in,dnedr_in,te_in,dtedr_in,ni_in,dnidr_in,ti_in,dtidr_in) + USE pprof_pass + IMPLICIT NONE + REAL(rknd), INTENT(IN) :: ne_in + REAL(rknd), INTENT(IN) :: dnedr_in + REAL(rknd), INTENT(IN) :: te_in + REAL(rknd), INTENT(IN) :: dtedr_in + REAL(rknd), DIMENSION(num_ion_species), INTENT(IN) :: ni_in + REAL(rknd), DIMENSION(num_ion_species), INTENT(IN) :: dnidr_in + REAL(rknd), DIMENSION(num_ion_species), INTENT(IN) :: ti_in + REAL(rknd), DIMENSION(num_ion_species), INTENT(IN) :: dtidr_in + ne = ne_in + dnedr = dnedr_in + te = te_in + dtedr = dtedr_in + ni = ni_in + dnidr = dnidr_in + ti = ni_in + dtidr = dnidr_in + RETURN + END SUBROUTINE penta_set_pprof + + SUBROUTINE penta_set_DKES_star(nc,ne,cmul_in,efield_in,D11_in,D13_in,D33_in) + USE vmec_var_pass, ONLY: Bsq + USE coeff_var_pass + IMPLICIT NONE + INTEGER(iknd), INTENT(in) :: nc,ne + REAL(rknd), DIMENSION(nc), INTENT(IN) :: cmul_in + REAL(rknd), DIMENSION(ne), INTENT(IN) :: efield_in + REAL(rknd), DIMENSION(nc,ne), INTENT(IN) :: D11_in + REAL(rknd), DIMENSION(nc,ne), INTENT(IN) :: D13_in + REAL(rknd), DIMENSION(nc,ne), INTENT(IN) :: D33_in + INTEGER :: je, ic + REAL(rknd) :: D33_Spitzer + num_c = nc + num_e = ne + IF (ALLOCATED(cmul)) DEALLOCATE(cmul) + IF (ALLOCATED(efield)) DEALLOCATE(efield) + IF (ALLOCATED(D11_mat)) DEALLOCATE(D11_mat) + IF (ALLOCATED(D13_mat)) DEALLOCATE(D13_mat) + IF (ALLOCATED(D31_mat)) DEALLOCATE(D31_mat) + IF (ALLOCATED(D33_mat)) DEALLOCATE(D33_mat) + ALLOCATE(cmul(num_c)) + ALLOCATE(efield(num_e)) + ALLOCATE(D11_mat(num_c,num_e)) + ALLOCATE(D13_mat(num_c,num_e)) + ALLOCATE(D31_mat(num_c,num_e)) + ALLOCATE(D33_mat(num_c,num_e)) + cmul = cmul_in + efield = efield_in + D11_mat = D11_in + D13_mat = D13_in + D33_mat = D33_in + ! Create D31 using Onsager symmetry + D31_mat = -D13_mat + IF (Add_Spitzer_to_D33) THEN + DO je = 1, num_e + DO ic = 1, num_c + ! Calculate D33* Spitzer = (2/3)*/cmul + D33_Spitzer = (2._rknd/3._rknd)*Bsq/cmul(ic) + ! Calculate D33*(Physical)=D33*(Spitzer)-D33* + D33_mat(ic,je) = D33_Spitzer - D33_mat(ic,je) + ENDDO + ENDDO + END IF + END SUBROUTINE penta_set_DKES_star + + SUBROUTINE penta_set_beam(beam_force_in) + Use pprof_pass, ONLY: beam_force + IMPLICIT NONE + REAL(rknd), INTENT(IN) :: beam_force_in + beam_force = beam_force_in + RETURN + END SUBROUTINE penta_set_beam + + SUBROUTINE penta_set_U2(U2_in) + USE coeff_var_pass, ONLY: D11_mat,cmul,num_c + IMPLICIT NONE + REAL(rknd), INTENT(IN), OPTIONAL :: U2_in + U2=1.5d0*D11_mat(num_c,1)/cmul(num_c) + IF (PRESENT(U2_in)) U2 = U2_in + RETURN + END SUBROUTINE penta_set_U2 + SUBROUTINE penta_fit_DXX_coef USE coeff_var_pass USE PENTA_subroutines, ONLY : fit_coeffs @@ -723,7 +814,7 @@ SUBROUTINE penta_run_1_init CALL read_penta_run_params_namelist('run_params',istat) CALL penta_init_commandline CALL penta_allocate_species - CALL penta_read_input_files + CALL penta_read_input_files(.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.) CALL penta_screen_info CALL penta_allocate_dkescoeff CALL penta_fit_DXX_coef From 5636ea4f84ceb104156eea8f24e6cfe56c49e7a2 Mon Sep 17 00:00:00 2001 From: lazersos Date: Fri, 18 Oct 2024 15:43:56 +0200 Subject: [PATCH 28/47] PENTA: Fixed a variable name to reflect dVds value not Volume. --- PENTA/Sources/penta_interface_mod.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/PENTA/Sources/penta_interface_mod.f90 b/PENTA/Sources/penta_interface_mod.f90 index 702d40be5..8cee618db 100644 --- a/PENTA/Sources/penta_interface_mod.f90 +++ b/PENTA/Sources/penta_interface_mod.f90 @@ -458,13 +458,13 @@ SUBROUTINE penta_read_input_files(lvmec,lpprof,ldkes,lbeam,lU2) RETURN END SUBROUTINE penta_read_input_files - SUBROUTINE penta_set_eq_data(rho,A_in,R_in,V_in,chip_in,phip_in,iota_in,bth_in,bze_in,Bsq_in) + SUBROUTINE penta_set_eq_data(rho,A_in,R_in,Vp_in,chip_in,phip_in,iota_in,bth_in,bze_in,Bsq_in) USE vmec_var_pass IMPLICIT NONE REAL(rknd), INTENT(IN) :: rho REAL(rknd), INTENT(IN) :: A_in REAL(rknd), INTENT(IN) :: R_in - REAL(rknd), INTENT(IN) :: V_in + REAL(rknd), INTENT(IN) :: Vp_in REAL(rknd), INTENT(IN) :: chip_in REAL(rknd), INTENT(IN) :: phip_in REAL(rknd), INTENT(IN) :: iota_in @@ -483,7 +483,7 @@ SUBROUTINE penta_set_eq_data(rho,A_in,R_in,V_in,chip_in,phip_in,iota_in,bth_in,b btheta = bth_in bzeta = bze_in Bsq = Bsq_in - vol_p = V_in + vol_p = Vp_in ! Note that vp from VMEC comes normalized by 4pi^2 ! Therefore we need to denormalize it vol_p = TWOPI*TWOPI*vol_p From b42c44d2d3c6aa9edea5176e9e3a8bcd08c71ef1 Mon Sep 17 00:00:00 2001 From: lazersos Date: Fri, 18 Oct 2024 15:44:45 +0200 Subject: [PATCH 29/47] THRIFT: Added buco, bvco, and bsq splines to support PENTA integration. --- THRIFT/Sources/thrift_equil.f90 | 3 ++- THRIFT/Sources/thrift_load_vmec.f90 | 30 +++++++++++++++++++++++++++++ 2 files changed, 32 insertions(+), 1 deletion(-) diff --git a/THRIFT/Sources/thrift_equil.f90 b/THRIFT/Sources/thrift_equil.f90 index e399c10aa..94bc2a236 100644 --- a/THRIFT/Sources/thrift_equil.f90 +++ b/THRIFT/Sources/thrift_equil.f90 @@ -27,7 +27,8 @@ MODULE thrift_equil ! Spline helpers INTEGER :: bcs1(2) - TYPE(EZspline1_r8) :: iota_spl, phip_spl, vp_spl + TYPE(EZspline1_r8) :: iota_spl, phip_spl, vp_spl, bu_spl, bv_spl, & + bsq_spl END MODULE thrift_equil diff --git a/THRIFT/Sources/thrift_load_vmec.f90 b/THRIFT/Sources/thrift_load_vmec.f90 index c6b50bf36..d1ed29d88 100644 --- a/THRIFT/Sources/thrift_load_vmec.f90 +++ b/THRIFT/Sources/thrift_load_vmec.f90 @@ -230,6 +230,36 @@ SUBROUTINE thrift_load_vmec CALL EZspline_setup(vp_spl,vp,iflag,EXACT_DIM=.true.) IF (iflag /=0) CALL handle_err(EZSPLINE_ERR,'thrift_load_vmec: vp',iflag) + ! BU + bcs1=(/ 0, 0/) + IF (EZspline_allocated(bu_spl)) CALL EZspline_free(bu_spl,iflag) + CALL EZspline_init(bu_spl,ns,bcs1,iflag) + IF (iflag /=0) CALL handle_err(EZSPLINE_ERR,'thrift_load_vmec: bu',iflag) + bu_spl%isHermite = 0 + FORALL (k=1:ns) bu_spl%x1(k) = sqrt(DBLE(k-1)/DBLE(ns-1)) + CALL EZspline_setup(bu_spl,buco,iflag,EXACT_DIM=.true.) + IF (iflag /=0) CALL handle_err(EZSPLINE_ERR,'thrift_load_vmec: bu',iflag) + + ! BV + bcs1=(/ 0, 0/) + IF (EZspline_allocated(bv_spl)) CALL EZspline_free(bv_spl,iflag) + CALL EZspline_init(bv_spl,ns,bcs1,iflag) + IF (iflag /=0) CALL handle_err(EZSPLINE_ERR,'thrift_load_vmec: bv',iflag) + bv_spl%isHermite = 0 + FORALL (k=1:ns) bv_spl%x1(k) = sqrt(DBLE(k-1)/DBLE(ns-1)) + CALL EZspline_setup(bv_spl,bvco,iflag,EXACT_DIM=.true.) + IF (iflag /=0) CALL handle_err(EZSPLINE_ERR,'thrift_load_vmec: bv',iflag) + + ! BSQ + bcs1=(/ 0, 0/) + IF (EZspline_allocated(bsq_spl)) CALL EZspline_free(bsq_spl,iflag) + CALL EZspline_init(bsq_spl,ns,bcs1,iflag) + IF (iflag /=0) CALL handle_err(EZSPLINE_ERR,'thrift_load_vmec: bsq',iflag) + bsq_spl%isHermite = 0 + FORALL (k=1:ns) bsq_spl%x1(k) = sqrt(DBLE(k-1)/DBLE(ns-1)) + CALL EZspline_setup(bsq_spl,bdotb,iflag,EXACT_DIM=.true.) + IF (iflag /=0) CALL handle_err(EZSPLINE_ERR,'thrift_load_vmec: bsq',iflag) + RETURN !---------------------------------------------------------------------- From 9ba34b7a36a110d2c6a162b9b15d0cd36ab6028b Mon Sep 17 00:00:00 2001 From: lazersos Date: Fri, 18 Oct 2024 15:45:07 +0200 Subject: [PATCH 30/47] THRIFT: PENTA interface updated. --- THRIFT/Sources/thrift_penta.f90 | 63 ++++++++++++++++++++++++++------- 1 file changed, 50 insertions(+), 13 deletions(-) diff --git a/THRIFT/Sources/thrift_penta.f90 b/THRIFT/Sources/thrift_penta.f90 index e63620b83..4fa6ab5cd 100644 --- a/THRIFT/Sources/thrift_penta.f90 +++ b/THRIFT/Sources/thrift_penta.f90 @@ -28,8 +28,11 @@ SUBROUTINE thrift_penta(lscreen,iflag) !----------------------------------------------------------------------- ! Local Variables !----------------------------------------------------------------------- - INTEGER :: ns_dkes, k - REAL(rprec) :: s + INTEGER :: ns_dkes, k, ier, j, i, ncstar, nestar + REAL(rprec) :: s, rho, iota, phip, chip, btheta, bzeta, bsq, vp, & + te, ne, dtedrho, dnedrho + REAL(rprec), DIMENSION(num_ion_species) :: ni,ti, dtidrho, dnidrho + REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: D11, D13, D33 !----------------------------------------------------------------------- ! BEGIN SUBROUTINE !----------------------------------------------------------------------- @@ -46,18 +49,52 @@ SUBROUTINE thrift_penta(lscreen,iflag) ! Break up work CALL MPI_CALC_MYRANGE(MPI_COMM_MYWORLD,1,ns_dkes,mystart,myend) DO k = mystart,myend - ! This part mimics penta_run_1 but without reading the namelists - Er_min = -250.0 - Er_max = 250.0 - js = DKES_K(k) - i_append = 1 - B_Eprl = 0.0 - Smax = 1 - coeff_ext = '' - run_ident = '' - pprof_char = '' + ! Calc some information + s = DKES_K(k)/ns_equil ! <- Need to see if we have this number + rho = sqrt(s) + ier = 0 + CALL EZSpline_interp(iota_spl, rho, iota, ier) + ier = 0 + CALL EZSpline_interp(phip_spl, rho, phip, ier) + ! PENTA wants d/ds note that this won't work if rho=0 + phip = 0.5*phip/rho + chip = iota * phip + ier = 0 + CALL EZSpline_interp(bu_spl, rho, btheta, ier) + ier = 0 + CALL EZSpline_interp(bv_spl, rho, bzeta, ier) + ier = 0 + CALL EZSpline_interp(bsq_spl, rho, bsq, ier) + ier = 0 + CALL EZSpline_interp(vp_spl, rho, vp, ier) + ! Vp = dV/dPHI need dVds with VMEC normalization + vp = vp * phip /(pi2*pi2) + ! Profiles + CALL get_prof_te(rho, THRIFT_T(mytimestep), te) + CALL get_prof_ne(rho, THRIFT_T(mytimestep), ne) + CALL get_prof_teprime(rho, THRIFT_T(mytimestep), dtedrho) + CALL get_prof_neprime(rho, THRIFT_T(mytimestep), dnedrho) + DO j = 1, num_ion_species + CALL get_prof_ti(rho, THRIFT_T(mytimestep), j, ti(j)) + CALL get_prof_ni(rho, THRIFT_T(mytimestep), j, ni(j)) + CALL get_prof_tiprime(rho, THRIFT_T(mytimestep), j, dtidrho(j)) + CALL get_prof_niprime(rho, THRIFT_T(mytimestep), j, dnidrho(j)) + ! DKES Data + ncstar = COUNT(DKES_NUSTAR < 1E10) + nestar = COUNT(DKES_ERSTAR < 1E10) + + ! PENTA + CALL PENTA_SET_ION_PARAMS(nion_prof, Zatom_prof, mass_prof/AMU) + CALL PENTA_SET_COMMANDLINE(-250.0,250.0,DKES_K(k),1,0.0,1,'','','') CALL PENTA_ALLOCATE_SPECIES - CALL PENTA_READ_INPUT_FILES + CALL PENTA_SET_EQ_DATA(rho,eq_Aminor,eq_Rmajor,vp,chip,phip,iota,btheta,bzeta,bsq) + CALL PENTA_SET_PPROF(ne,dnedrho,te,dtedrho,ni,dnidrho,ti,dtidrho) + !SUBROUTINE penta_set_DKES_star(nc,ne,cmul_in,efield_in,D11_in,D13_in,D33_in) + CALL PENTA_SET_DKES_STAR(ncstar,nestar,DKES_NUSTAR(1:ncstar),DKES_ERSTAR(1:nestar), & + D11,D13,D33) + CALL PENTA_SET_BEAM(0.0) ! Zero becasue we don't read + CALL PENTA_SET_U2() ! Leave blank for default value + CALL PENTA_READ_INPUT_FILES(.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.) CALL PENTA_SCREEN_INFO CALL PENTA_ALLOCATE_DKESCOEFF CALL PENTA_FIT_DXX_COEF From ba14122f043d6de81480969a2d7a42f39b3667b7 Mon Sep 17 00:00:00 2001 From: lazersos Date: Fri, 18 Oct 2024 15:58:27 +0200 Subject: [PATCH 31/47] PENTA: Fixed type of Smax in interface routine. --- PENTA/Sources/penta_interface_mod.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PENTA/Sources/penta_interface_mod.f90 b/PENTA/Sources/penta_interface_mod.f90 index 8cee618db..fd12a703c 100644 --- a/PENTA/Sources/penta_interface_mod.f90 +++ b/PENTA/Sources/penta_interface_mod.f90 @@ -282,7 +282,7 @@ SUBROUTINE penta_set_commandline(Er_min_in,Er_max_in,js_in,i_append_in,B_Eprl_in INTEGER(iknd), INTENT(IN) :: js_in INTEGER(iknd), INTENT(IN) :: i_append_in REAL(rknd), INTENT(IN) :: B_Eprl_in - REAL(rknd), INTENT(IN) :: Smax_in + INTEGER(iknd), INTENT(IN) :: Smax_in CHARACTER(LEN=*), INTENT(IN) :: ext_in CHARACTER(LEN=*), INTENT(IN) :: run_in CHARACTER(LEN=*), INTENT(IN) :: pprof_in From 32d1db3faacea26efd96737399beda979fff23a0 Mon Sep 17 00:00:00 2001 From: lazersos Date: Fri, 18 Oct 2024 15:59:03 +0200 Subject: [PATCH 32/47] THRIFT: Added VMEC ns helper --- THRIFT/Sources/thrift_equil.f90 | 2 +- THRIFT/Sources/thrift_load_vmec.f90 | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/THRIFT/Sources/thrift_equil.f90 b/THRIFT/Sources/thrift_equil.f90 index 94bc2a236..41147cfc1 100644 --- a/THRIFT/Sources/thrift_equil.f90 +++ b/THRIFT/Sources/thrift_equil.f90 @@ -23,7 +23,7 @@ MODULE thrift_equil IMPLICIT NONE REAL(rprec) :: bt0, el0, iota0, eq_beta, eq_Aminor, eq_Rmajor, & - eq_phiedge, eq_volume + eq_phiedge, eq_volume, ns_eq ! Spline helpers INTEGER :: bcs1(2) diff --git a/THRIFT/Sources/thrift_load_vmec.f90 b/THRIFT/Sources/thrift_load_vmec.f90 index d1ed29d88..184cf42ca 100644 --- a/THRIFT/Sources/thrift_load_vmec.f90 +++ b/THRIFT/Sources/thrift_load_vmec.f90 @@ -184,6 +184,7 @@ SUBROUTINE thrift_load_vmec ! Helpers b0 = RBtor0/Rmajor iota0 = iotaf(1) + ns_eq = ns ! We want dV/dPhi = dV/ds*ds/dPhi = dV/ds / (dPhi/ds) ! Note that Vp is missing a factor of 4*pi*pi From 15c86acc12d64d67c0f221ccc4bbdda7cef9d703 Mon Sep 17 00:00:00 2001 From: lazersos Date: Fri, 18 Oct 2024 15:59:40 +0200 Subject: [PATCH 33/47] THRIFT: Made AMU available through the profiles module. --- THRIFT/Sources/thrift_profiles_mod.f90 | 2 ++ THRIFT/Sources/thrift_run_nbcd.f90 | 1 - 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/THRIFT/Sources/thrift_profiles_mod.f90 b/THRIFT/Sources/thrift_profiles_mod.f90 index 3bf2b5132..23f21bea8 100644 --- a/THRIFT/Sources/thrift_profiles_mod.f90 +++ b/THRIFT/Sources/thrift_profiles_mod.f90 @@ -26,6 +26,8 @@ MODULE thrift_profiles_mod INTEGER :: win_raxis_prof, win_taxis_prof, win_NE3D, win_TE3D, & win_NI4D, win_TI4D, win_hr, win_ht, win_hri, win_hti, & win_Matom_prof, win_Zatom_prof, win_P3D + + REAL(rprec), PARAMETER :: AMU = 1.66053906892D-27 !----------------------------------------------------------------------- ! Input Namelists ! NONE diff --git a/THRIFT/Sources/thrift_run_nbcd.f90 b/THRIFT/Sources/thrift_run_nbcd.f90 index e0db38e0c..1595fd749 100644 --- a/THRIFT/Sources/thrift_run_nbcd.f90 +++ b/THRIFT/Sources/thrift_run_nbcd.f90 @@ -17,7 +17,6 @@ SUBROUTINE thrift_run_nbcd ! i Loop index !----------------------------------------------------------------------- IMPLICIT NONE - REAL(rprec), PARAMETER :: amu = 1.66053906660D-27 INTEGER :: i REAL(rprec) :: mass, NBI_A, NBI_E, NBI_Z, NBI_alpha, NBI_theta, NBI_Pdense ! probably should be inputs REAL(rprec) :: ne, Te, Ti, Zeff, InvAspect, G, NBI_V, rho, timenow, jopb From 8ceb5cef1eb3cc2a05b4c222747794db6a7041ae Mon Sep 17 00:00:00 2001 From: lazersos Date: Fri, 18 Oct 2024 16:00:42 +0200 Subject: [PATCH 34/47] THRIFT: PENTA added, missing DXX coefficient handling. --- THRIFT/ObjectList | 1 + THRIFT/Sources/thrift_penta.f90 | 13 ++++++------- THRIFT/THRIFT.dep | 10 ++++++++++ 3 files changed, 17 insertions(+), 7 deletions(-) diff --git a/THRIFT/ObjectList b/THRIFT/ObjectList index 1da248bf7..f237bd14b 100644 --- a/THRIFT/ObjectList +++ b/THRIFT/ObjectList @@ -1,4 +1,5 @@ ObjectFiles = \ +thrift_penta.o \ thrift_dkes.o \ thrift_diagno.o \ thrift_run_diagnostics.o \ diff --git a/THRIFT/Sources/thrift_penta.f90 b/THRIFT/Sources/thrift_penta.f90 index 4fa6ab5cd..26d972821 100644 --- a/THRIFT/Sources/thrift_penta.f90 +++ b/THRIFT/Sources/thrift_penta.f90 @@ -9,7 +9,6 @@ SUBROUTINE thrift_penta(lscreen,iflag) ! Libraries !----------------------------------------------------------------------- USE thrift_runtime - USE thrift_input_mod USE thrift_vars, nrho_thrift => nrho USE thrift_profiles_mod USE thrift_equil @@ -28,7 +27,7 @@ SUBROUTINE thrift_penta(lscreen,iflag) !----------------------------------------------------------------------- ! Local Variables !----------------------------------------------------------------------- - INTEGER :: ns_dkes, k, ier, j, i, ncstar, nestar + INTEGER :: ns_dkes, k, ier, j, i, ncstar, nestar, mystart, myend REAL(rprec) :: s, rho, iota, phip, chip, btheta, bzeta, bsq, vp, & te, ne, dtedrho, dnedrho REAL(rprec), DIMENSION(num_ion_species) :: ni,ti, dtidrho, dnidrho @@ -36,7 +35,6 @@ SUBROUTINE thrift_penta(lscreen,iflag) !----------------------------------------------------------------------- ! BEGIN SUBROUTINE !----------------------------------------------------------------------- - IMPLICIT NONE IF (iflag < 0) RETURN IF (lscreen) WRITE(6,'(a)') ' -------------------- NEOCLASSICAL FLUX USING PENTA -------------------' IF (lvmec) THEN @@ -50,7 +48,7 @@ SUBROUTINE thrift_penta(lscreen,iflag) CALL MPI_CALC_MYRANGE(MPI_COMM_MYWORLD,1,ns_dkes,mystart,myend) DO k = mystart,myend ! Calc some information - s = DKES_K(k)/ns_equil ! <- Need to see if we have this number + s = DKES_K(k)/ns_eq rho = sqrt(s) ier = 0 CALL EZSpline_interp(iota_spl, rho, iota, ier) @@ -79,20 +77,21 @@ SUBROUTINE thrift_penta(lscreen,iflag) CALL get_prof_ni(rho, THRIFT_T(mytimestep), j, ni(j)) CALL get_prof_tiprime(rho, THRIFT_T(mytimestep), j, dtidrho(j)) CALL get_prof_niprime(rho, THRIFT_T(mytimestep), j, dnidrho(j)) + END DO ! DKES Data ncstar = COUNT(DKES_NUSTAR < 1E10) nestar = COUNT(DKES_ERSTAR < 1E10) ! PENTA - CALL PENTA_SET_ION_PARAMS(nion_prof, Zatom_prof, mass_prof/AMU) - CALL PENTA_SET_COMMANDLINE(-250.0,250.0,DKES_K(k),1,0.0,1,'','','') + CALL PENTA_SET_ION_PARAMS(nion_prof, DBLE(Zatom_prof), Matom_prof/AMU) + CALL PENTA_SET_COMMANDLINE(-250.0_rprec,250.0_rprec,DKES_K(k),1,0.0_rprec,1,'','','') CALL PENTA_ALLOCATE_SPECIES CALL PENTA_SET_EQ_DATA(rho,eq_Aminor,eq_Rmajor,vp,chip,phip,iota,btheta,bzeta,bsq) CALL PENTA_SET_PPROF(ne,dnedrho,te,dtedrho,ni,dnidrho,ti,dtidrho) !SUBROUTINE penta_set_DKES_star(nc,ne,cmul_in,efield_in,D11_in,D13_in,D33_in) CALL PENTA_SET_DKES_STAR(ncstar,nestar,DKES_NUSTAR(1:ncstar),DKES_ERSTAR(1:nestar), & D11,D13,D33) - CALL PENTA_SET_BEAM(0.0) ! Zero becasue we don't read + CALL PENTA_SET_BEAM(0.0_rprec) ! Zero becasue we don't read CALL PENTA_SET_U2() ! Leave blank for default value CALL PENTA_READ_INPUT_FILES(.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.) CALL PENTA_SCREEN_INFO diff --git a/THRIFT/THRIFT.dep b/THRIFT/THRIFT.dep index 018053f37..53398edc0 100644 --- a/THRIFT/THRIFT.dep +++ b/THRIFT/THRIFT.dep @@ -135,6 +135,16 @@ thrift_paraexe.o : \ $(VMEC_DIR)/$(LOCTYPE)/parallel_vmec_module.o \ $(VMEC_DIR)/$(LOCTYPE)/vmec_params.o +thrift_penta.o : \ + thrift_runtime.o \ + thrift_vars.o \ + thrift_profiles_mod.o \ + thrift_equil.o \ + thrift_funcs.o \ + $(LIB_DIR)/$(LOCTYPE)/mpi_inc.o \ + $(LIB_DIR)/$(LOCTYPE)/mpi_params.o \ + $(PENTA_DIR)/$(LOCTYPE)/penta_interface_mod.o + thrift_profiles_mod.o : \ thrift_runtime.o \ $(LIB_DIR)/$(LOCTYPE)/stel_kinds.o \ From abfab33224eca67b003e61e142423be8e9598b8b Mon Sep 17 00:00:00 2001 From: lazersos Date: Fri, 18 Oct 2024 17:21:32 +0200 Subject: [PATCH 35/47] DKES: Added some helper vars for THRIFT. --- DKES/Sources/General/dkes_realspace.f | 1 + 1 file changed, 1 insertion(+) diff --git a/DKES/Sources/General/dkes_realspace.f b/DKES/Sources/General/dkes_realspace.f index 362ccdd2b..0d297d6d2 100644 --- a/DKES/Sources/General/dkes_realspace.f +++ b/DKES/Sources/General/dkes_realspace.f @@ -27,6 +27,7 @@ MODULE dkes_realspace 1 DKES_L11p, DKES_L33p, DKES_L31p, 2 DKES_L11m, DKES_L33m, DKES_L31m, 3 DKES_scal11, DKES_scal33, DKES_scal31 + INTEGER :: DKES_NK, DKES_NC, DKES_NE CONTAINS From 1ea2341d3b23a21e36d167ae31e418cb94cadbb0 Mon Sep 17 00:00:00 2001 From: lazersos Date: Fri, 18 Oct 2024 17:22:41 +0200 Subject: [PATCH 36/47] THRIFT: DKES and PENTA now integrated. --- THRIFT/Sources/thrift_boozer.f90 | 10 +- THRIFT/Sources/thrift_dkes.f90 | 22 +++++ THRIFT/Sources/thrift_equil.f90 | 4 +- THRIFT/Sources/thrift_init.f90 | 16 +++- THRIFT/Sources/thrift_paraexe.f90 | 10 ++ THRIFT/Sources/thrift_penta.f90 | 10 +- THRIFT/Sources/thrift_run_bootstrap.f90 | 121 ++++++++++++++---------- THRIFT/Sources/thrift_runtime.f90 | 4 + THRIFT/THRIFT.dep | 5 +- 9 files changed, 135 insertions(+), 67 deletions(-) diff --git a/THRIFT/Sources/thrift_boozer.f90 b/THRIFT/Sources/thrift_boozer.f90 index 896f24e00..4921c7e3d 100644 --- a/THRIFT/Sources/thrift_boozer.f90 +++ b/THRIFT/Sources/thrift_boozer.f90 @@ -63,11 +63,11 @@ SUBROUTINE thrift_boozer(lscreen,iflag) CALL MPI_BCAST(nboz_xboozer,1,MPI_INTEGER,master,MPI_COMM_MYWORLD,ierr_mpi) CALL MPI_BCAST(ns_vmec,1,MPI_INTEGER,master,MPI_COMM_MYWORLD,ierr_mpi) #endif - ! We need to setup the lsurf_boz array - IF (ALLOCATED(lsurf_boz)) DEALLOCATE(lsurf_boz) - ALLOCATE(lsurf_boz(ns_vmec)) - lsurf_boz = .FALSE. - lsurf_boz(2:ns_vmec) = .TRUE. ! We can probaly do less surfaces in the future + ! We need to setup the lsurf_boz array (this in now calculated outside this routine) + !IF (ALLOCATED(lsurf_boz)) DEALLOCATE(lsurf_boz) + !ALLOCATE(lsurf_boz(ns_vmec)) + !lsurf_boz = .FALSE. + !lsurf_boz(2:ns_vmec) = .TRUE. ! We can probaly do less surfaces in the future num_booz = -COUNT(lsurf_boz) ! This tricks the code into not reading the wout or input file ! Now read the wout file ier = 0 diff --git a/THRIFT/Sources/thrift_dkes.f90 b/THRIFT/Sources/thrift_dkes.f90 index f55a45fb4..e3a38ebda 100644 --- a/THRIFT/Sources/thrift_dkes.f90 +++ b/THRIFT/Sources/thrift_dkes.f90 @@ -17,6 +17,7 @@ SUBROUTINE thrift_dkes(lscreen,iflag) USE thrift_funcs USE mpi_params USE mpi_inc + USE mpi_sharmem USE read_wout_mod, ONLY: read_wout_file, read_wout_deallocate USE read_boozer_mod, ONLY: bcast_boozer_vars ! DKES LIBRARIES @@ -84,6 +85,9 @@ SUBROUTINE thrift_dkes(lscreen,iflag) DKES_scal11=0.0; DKES_scal33=0.0; DKES_scal31=0.0 l = 1 ! Allocate Helper arrays to loop over dkes runs + DKES_NK = COUNT(DKES_K>0) + DKES_NC = COUNT(DKES_Nustar<1E10) + DKES_NE = COUNT(DKES_Erstar<1E10) ALLOCATE(ik_dkes(nruns_dkes),Earr_dkes(nruns_dkes),nuarr_dkes(nruns_dkes)) IF (myworkid == master) THEN DO k = 1, DKES_NS_MAX @@ -288,6 +292,24 @@ SUBROUTINE thrift_dkes(lscreen,iflag) IF (lscreen) WRITE(6,'(a)') ' ----------------- DKES CALCULATION (DONE) ----------------' IF (myworkid .ne. master) CALL read_wout_deallocate + ! Sort out the DKES runs for later use + IF (ASSOCIATED(DKES_D11)) CALL mpidealloc(DKES_D11,win_dkes_d11) + IF (ASSOCIATED(DKES_D31)) CALL mpidealloc(DKES_D31,win_dkes_d31) + IF (ASSOCIATED(DKES_D33)) CALL mpidealloc(DKES_D33,win_dkes_d33) + CALL mpialloc(DKES_D11, DKES_NK, DKES_NC, DKES_NE, myid_sharmem, 0, MPI_COMM_MYWORLD, win_dkes_d11) + CALL mpialloc(DKES_D31, DKES_NK, DKES_NC, DKES_NE, myid_sharmem, 0, MPI_COMM_MYWORLD, win_dkes_d31) + CALL mpialloc(DKES_D33, DKES_NK, DKES_NC, DKES_NE, myid_sharmem, 0, MPI_COMM_MYWORLD, win_dkes_d33) + IF (myworkid .eq. master) THEN + DO l = 1, nruns_dkes + i = MOD(l-1,DKES_NK)+1 + j = MOD(l-1,DKES_NK*DKES_NE) + j = FLOOR(REAL(j) / REAL(DKES_NK))+1 + k = CEILING(REAL(l) / REAL(DKES_NK*DKES_NE)) + DKES_D11(i,j,k) = (DKES_L11p(l) + DKES_L11m(l))*0.5 + DKES_D31(i,j,k) = (DKES_L31p(l) + DKES_L31m(l))*0.5 + DKES_D33(i,j,k) = (DKES_L33p(l) + DKES_L33m(l))*0.5 + END DO + END IF END IF IF (lscreen) WRITE(6,'(a)') ' ------------------- DKES NEOCLASSICAL CALCULATION DONE ---------------------' RETURN diff --git a/THRIFT/Sources/thrift_equil.f90 b/THRIFT/Sources/thrift_equil.f90 index 41147cfc1..950474824 100644 --- a/THRIFT/Sources/thrift_equil.f90 +++ b/THRIFT/Sources/thrift_equil.f90 @@ -21,9 +21,9 @@ MODULE thrift_equil ! vp_spl Spline of dV/dphi !------------------------------------------------------------------- IMPLICIT NONE - + INTEGER :: ns_eq REAL(rprec) :: bt0, el0, iota0, eq_beta, eq_Aminor, eq_Rmajor, & - eq_phiedge, eq_volume, ns_eq + eq_phiedge, eq_volume ! Spline helpers INTEGER :: bcs1(2) diff --git a/THRIFT/Sources/thrift_init.f90 b/THRIFT/Sources/thrift_init.f90 index 99d6dad86..70032608e 100644 --- a/THRIFT/Sources/thrift_init.f90 +++ b/THRIFT/Sources/thrift_init.f90 @@ -13,7 +13,9 @@ SUBROUTINE thrift_init USE thrift_input_mod USE thrift_vars USE thrift_profiles_mod - USE diagno_input_mod, ONLY: read_diagno_input + USE diagno_input_mod, ONLY: read_diagno_input + USE penta_interface_mod, ONLY: init_penta_input, & + read_penta_run_params_namelist USE safe_open_mod USE mpi_params USE mpi_inc @@ -61,8 +63,6 @@ SUBROUTINE thrift_init WRITE(6,'(A11,F8.4)') ' TEND: ', tend END IF - ! Create the worker pool - ! Grid allocations CALL mpialloc(THRIFT_RHO, nrho, myid_sharmem, 0, MPI_COMM_SHARMEM, win_thrift_rho) CALL mpialloc(THRIFT_RHOFULL, nrho+2, myid_sharmem, 0, MPI_COMM_SHARMEM, win_thrift_rhofull) @@ -142,6 +142,16 @@ SUBROUTINE thrift_init STOP END IF CLOSE(iunit) + CASE('dkespenta') + ier = 0 + CALL init_penta_input + CALL read_penta_run_params_namelist('input.'//TRIM(id_string),ier) + IF (ier < 0 .and. myid == master) THEN + WRITE(6,*) '!!!!!!!!!!!!ERRROR!!!!!!!!!!!!!!' + WRITE(6,*) ' RUN_PARAMS Namelist not found ' + WRITE(6,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + STOP + END IF CASE('sfincs') END SELECT diff --git a/THRIFT/Sources/thrift_paraexe.f90 b/THRIFT/Sources/thrift_paraexe.f90 index 986a2d9f5..0275f7871 100644 --- a/THRIFT/Sources/thrift_paraexe.f90 +++ b/THRIFT/Sources/thrift_paraexe.f90 @@ -185,6 +185,16 @@ SUBROUTINE thrift_paraexe(in_parameter_1,in_parameter_2,lscreen_local) ier = 0 CALL thrift_diagno(lscreen_local,ier) ier_paraexe = ier + CASE('dkes') + proc_string = file_str + ier = 0 + CALL thrift_dkes(lscreen_local,ier) + ier_paraexe = ier + CASE('penta') + proc_string = file_str + ier = 0 + CALL thrift_penta(lscreen_local,ier) + ier_paraexe = ier CASE('sfincs') ! proc_string = file_str ! ier = 0 diff --git a/THRIFT/Sources/thrift_penta.f90 b/THRIFT/Sources/thrift_penta.f90 index 26d972821..968fff8d2 100644 --- a/THRIFT/Sources/thrift_penta.f90 +++ b/THRIFT/Sources/thrift_penta.f90 @@ -27,7 +27,8 @@ SUBROUTINE thrift_penta(lscreen,iflag) !----------------------------------------------------------------------- ! Local Variables !----------------------------------------------------------------------- - INTEGER :: ns_dkes, k, ier, j, i, ncstar, nestar, mystart, myend + INTEGER :: ns_dkes, k, ier, j, i, ncstar, nestar, mystart, myend, & + mysurf REAL(rprec) :: s, rho, iota, phip, chip, btheta, bzeta, bsq, vp, & te, ne, dtedrho, dnedrho REAL(rprec), DIMENSION(num_ion_species) :: ni,ti, dtidrho, dnidrho @@ -48,7 +49,8 @@ SUBROUTINE thrift_penta(lscreen,iflag) CALL MPI_CALC_MYRANGE(MPI_COMM_MYWORLD,1,ns_dkes,mystart,myend) DO k = mystart,myend ! Calc some information - s = DKES_K(k)/ns_eq + mysurf = DKES_K(k) + s = DBLE(mysurf)/DBLE(ns_eq) rho = sqrt(s) ier = 0 CALL EZSpline_interp(iota_spl, rho, iota, ier) @@ -88,9 +90,9 @@ SUBROUTINE thrift_penta(lscreen,iflag) CALL PENTA_ALLOCATE_SPECIES CALL PENTA_SET_EQ_DATA(rho,eq_Aminor,eq_Rmajor,vp,chip,phip,iota,btheta,bzeta,bsq) CALL PENTA_SET_PPROF(ne,dnedrho,te,dtedrho,ni,dnidrho,ti,dtidrho) - !SUBROUTINE penta_set_DKES_star(nc,ne,cmul_in,efield_in,D11_in,D13_in,D33_in) + ! Check here we have D31 from DKES but need D13 in PENTA CALL PENTA_SET_DKES_STAR(ncstar,nestar,DKES_NUSTAR(1:ncstar),DKES_ERSTAR(1:nestar), & - D11,D13,D33) + DKES_D11(mysurf,:,:),DKES_D31(mysurf,:,:),DKES_D33(mysurf,:,:)) CALL PENTA_SET_BEAM(0.0_rprec) ! Zero becasue we don't read CALL PENTA_SET_U2() ! Leave blank for default value CALL PENTA_READ_INPUT_FILES(.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.) diff --git a/THRIFT/Sources/thrift_run_bootstrap.f90 b/THRIFT/Sources/thrift_run_bootstrap.f90 index 8468e08b2..2c9e2bd4b 100644 --- a/THRIFT/Sources/thrift_run_bootstrap.f90 +++ b/THRIFT/Sources/thrift_run_bootstrap.f90 @@ -5,58 +5,75 @@ ! Description: This subroutine calls the relevant bootstrap ! current code. !----------------------------------------------------------------------- - SUBROUTINE thrift_run_bootstrap -!----------------------------------------------------------------------- -! Libraries -!----------------------------------------------------------------------- - USE thrift_runtime - USE thrift_vars - USE thrift_equil - USE thrift_profiles_mod - USE thrift_funcs -!----------------------------------------------------------------------- -! Local Variables -! ier Error flag -!----------------------------------------------------------------------- - IMPLICIT NONE - INTEGER :: i, ier - REAL(rprec) :: epsilon, pp1, pm1, ds - REAL(rprec), DIMENSION(:), ALLOCATABLE :: j_temp -!---------------------------------------------------------------------- -! BEGIN SUBROUTINE -!---------------------------------------------------------------------- - - ier = 0 - - ! Check to make sure we're not zero beta - IF (eq_beta == 0) THEN - THRIFT_JBOOT(:,mytimestep) = 0 - RETURN - END IF - - SELECT CASE(TRIM(bootstrap_type)) - CASE ('model','simple','test') - - ! j_BS = sqrt(epsilon) Rmajor *dp/dPhi - ! epsilon = a/R (inverse aspect ratio) - ! dp/dPhi : Pa/Wb (toroidal flux derivative) - ! dp/dPhi = dp/ds * ds/dPhi - ! = dp/ds / Phi_edge - ! j_BS = s*sqrt(epsilon)*Rmajor/Phi_edge*dp/ds - ! = s*sqrt(aminor*Rmajor)/Phi_edge*dp/ds - - epsilon = eq_Aminor/eq_Rmajor ! Inverse Aspect Ratio - THRIFT_JBOOT(:,mytimestep) = SQRT(eq_Aminor*eq_Rmajor)/eq_phiedge*THRIFT_S*THRIFT_PPRIME(:,mytimestep) - - CASE ('bootsj') - CALL thrift_paraexe('booz_xform',proc_string,lscreen_subcodes) - CALL thrift_paraexe('bootsj',proc_string,lscreen_subcodes) - CASE ('sfincs') - END SELECT +SUBROUTINE thrift_run_bootstrap + !----------------------------------------------------------------------- + ! Libraries + !----------------------------------------------------------------------- + USE thrift_runtime + USE thrift_vars + USE thrift_equil + USE thrift_profiles_mod + USE thrift_funcs + USE booz_params, ONLY: lsurf_boz + !----------------------------------------------------------------------- + ! Local Variables + ! ier Error flag + !----------------------------------------------------------------------- + IMPLICIT NONE + INTEGER :: i, ier + REAL(rprec) :: epsilon, pp1, pm1, ds + REAL(rprec), DIMENSION(:), ALLOCATABLE :: j_temp + !---------------------------------------------------------------------- + ! BEGIN SUBROUTINE + !---------------------------------------------------------------------- + + ier = 0 + ! Just always check it's deallocated + IF (ALLOCATED(lsurf_boz)) DEALLOCATE(lsurf_boz) + + ! Check to make sure we're not zero beta + IF (eq_beta == 0) THEN + THRIFT_JBOOT(:,mytimestep) = 0 RETURN -!---------------------------------------------------------------------- -! END SUBROUTINE -!---------------------------------------------------------------------- - END SUBROUTINE thrift_run_bootstrap + END IF + + SELECT CASE(TRIM(bootstrap_type)) + CASE ('model','simple','test') + + ! j_BS = sqrt(epsilon) Rmajor *dp/dPhi + ! epsilon = a/R (inverse aspect ratio) + ! dp/dPhi : Pa/Wb (toroidal flux derivative) + ! dp/dPhi = dp/ds * ds/dPhi + ! = dp/ds / Phi_edge + ! j_BS = s*sqrt(epsilon)*Rmajor/Phi_edge*dp/ds + ! = s*sqrt(aminor*Rmajor)/Phi_edge*dp/ds + + epsilon = eq_Aminor/eq_Rmajor ! Inverse Aspect Ratio + THRIFT_JBOOT(:,mytimestep) = SQRT(eq_Aminor*eq_Rmajor)/eq_phiedge*THRIFT_S*THRIFT_PPRIME(:,mytimestep) + + CASE ('bootsj') + ALLOCATE(lsurf_boz(ns_eq)) + lsurf_boz = .FALSE. + lsurf_boz(2:ns_eq) = .TRUE. + CALL thrift_paraexe('booz_xform',proc_string,lscreen_subcodes) + CALL thrift_paraexe('bootsj',proc_string,lscreen_subcodes) + CASE ('dkespenta') + ALLOCATE(lsurf_boz(ns_eq)) + lsurf_boz = .FALSE. + DO i = 1, DKES_NS_MAX + IF (DKES_K(i)<1) CYCLE + lsurf_boz(DKES_K(i)) = .TRUE. + END DO + CALL thrift_paraexe('booz_xform',proc_string,lscreen_subcodes) + CALL thrift_paraexe('dkes',proc_string,lscreen_subcodes) + CALL thrift_paraexe('penta',proc_string,lscreen_subcodes) + CASE ('sfincs') + END SELECT + + RETURN + !---------------------------------------------------------------------- + ! END SUBROUTINE + !---------------------------------------------------------------------- + END SUBROUTINE thrift_run_bootstrap diff --git a/THRIFT/Sources/thrift_runtime.f90 b/THRIFT/Sources/thrift_runtime.f90 index c6a4eab77..1c42290fe 100644 --- a/THRIFT/Sources/thrift_runtime.f90 +++ b/THRIFT/Sources/thrift_runtime.f90 @@ -78,6 +78,10 @@ MODULE thrift_runtime CHARACTER(256) :: id_string, prof_string, & nbcd_type, proc_string, magdiag_coil + LOGICAL, DIMENSION(MAXPROFLEN) :: lneed_boozer + REAL(rprec), DIMENSION(:,:,:), POINTER :: DKES_D11, DKES_D31, DKES_D33 + INTEGER :: win_dkes_d11, win_dkes_d31, win_dkes_d33 + REAL(rprec), PARAMETER :: THRIFT_VERSION = 0.50 !------------------------------------------------------------------- ! Subroutines diff --git a/THRIFT/THRIFT.dep b/THRIFT/THRIFT.dep index 53398edc0..2319d4038 100644 --- a/THRIFT/THRIFT.dep +++ b/THRIFT/THRIFT.dep @@ -78,6 +78,8 @@ thrift_init.o : \ thrift_runtime.o \ thrift_vars.o \ thrift_profiles_mod.o \ + $(LIB_DIR)/$(LOCTYPE)/diagno_input_mod.o \ + $(PENTA_DIR)/$(LOCTYPE)/penta_interface_mod.o \ $(LIB_DIR)/$(LOCTYPE)/mpi_inc.o \ $(LIB_DIR)/$(LOCTYPE)/mpi_params.o \ $(LIB_DIR)/$(LOCTYPE)/mpi_sharmem.o \ @@ -173,7 +175,8 @@ thrift_run_bootstrap.o : \ thrift_runtime.o \ thrift_equil.o \ thrift_vars.o \ - thrift_funcs.o + thrift_funcs.o \ + $(BOOZ_DIR)/$(LOCTYPE)/booz_params.o thrift_run_diagnostics.o : \ thrift_runtime.o \ From 1ac11b6c34dc9cd39f87e86112e75ec55b1d2d26 Mon Sep 17 00:00:00 2001 From: ajchcoelho Date: Wed, 13 Nov 2024 09:13:10 +0100 Subject: [PATCH 37/47] THRIFT: allow possibility of ntimesteps=1 (which corresponds to dt=0) --- THRIFT/Sources/thrift_init.f90 | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/THRIFT/Sources/thrift_init.f90 b/THRIFT/Sources/thrift_init.f90 index 8ea32d754..8c55662d8 100644 --- a/THRIFT/Sources/thrift_init.f90 +++ b/THRIFT/Sources/thrift_init.f90 @@ -233,7 +233,19 @@ SUBROUTINE thrift_init ENDIF ! Define grids - dt = (tend-tstart)/(ntimesteps-1) + IF( ntimesteps==1 ) THEN + dt = 0.0_rprec + ELSE IF( ntimesteps > 1) THEN + dt = (tend-tstart)/(ntimesteps-1) + ELSE + IF(lverb) THEN + WRITE(6,*) '!!!!!!!!!!!!ERRROR!!!!!!!!!!!!!!' + WRITE(6,*) ' ntimesteps < 1 ' + WRITE(6,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + STOP + END IF + END IF + IF (myid_sharmem == master) THEN FORALL(i = 1:nrho) THRIFT_RHO(i) = DBLE(i-0.5)/DBLE(nrho) ! (half) rho grid From 61826a5f9175cd9d158dde01f86d7d89650384b9 Mon Sep 17 00:00:00 2001 From: ajchcoelho Date: Wed, 13 Nov 2024 09:14:08 +0100 Subject: [PATCH 38/47] THRIFT: add comment when running Bootstrap --- THRIFT/Sources/thrift_evolve.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/THRIFT/Sources/thrift_evolve.f90 b/THRIFT/Sources/thrift_evolve.f90 index 8bd420ec7..763981d52 100644 --- a/THRIFT/Sources/thrift_evolve.f90 +++ b/THRIFT/Sources/thrift_evolve.f90 @@ -99,6 +99,7 @@ SUBROUTINE thrift_evolve CALL update_vars ! Calculate Bootstrap + IF (lverbj) WRITE(6,*) "Running Bootstrap code" CALL thrift_run_bootstrap THRIFT_JBOOT(:,mytimestep) = boot_factor*THRIFT_JBOOT(:,mytimestep) From 26aee6b3f418df0c22dd8219cd5f5c94430d3a7e Mon Sep 17 00:00:00 2001 From: ajchcoelho Date: Wed, 13 Nov 2024 09:28:17 +0100 Subject: [PATCH 39/47] THRIFT: printing of vprime was not being done correctly; not it is --- THRIFT/Sources/thrift_funcs.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/THRIFT/Sources/thrift_funcs.f90 b/THRIFT/Sources/thrift_funcs.f90 index a1df42a44..8dc9df92e 100644 --- a/THRIFT/Sources/thrift_funcs.f90 +++ b/THRIFT/Sources/thrift_funcs.f90 @@ -188,7 +188,7 @@ SUBROUTINE print_calc_vars() WRITE(6,*)' S VPRIME BAV BSQAV S11 ETAPARA PPRIME' WRITE(6,*)'' DO i = 1, nsj - WRITE(6,'(F5.3,3(1X,F7.3),3(1X,ES10.3))') & + WRITE(6,'(F5.3,1(1X,ES10.3),2(1X,F7.3),3(1X,ES10.3))') & THRIFT_S(i),THRIFT_VP(i,mytimestep),THRIFT_BAV(i,mytimestep),THRIFT_BSQAV(i,mytimestep), & THRIFT_S11(i,mytimestep),THRIFT_ETAPARA(i,mytimestep),THRIFT_PPRIME(i,mytimestep) END DO From ea70e71bf185e86cd622ebe19adf0e1bc73de5a6 Mon Sep 17 00:00:00 2001 From: ajchcoelho Date: Wed, 13 Nov 2024 15:25:19 +0100 Subject: [PATCH 40/47] DKES: remove prints --- DKES/Sources/General/printout.f | 2 -- 1 file changed, 2 deletions(-) diff --git a/DKES/Sources/General/printout.f b/DKES/Sources/General/printout.f index 9ed929a6f..2160ed7ac 100644 --- a/DKES/Sources/General/printout.f +++ b/DKES/Sources/General/printout.f @@ -145,7 +145,6 @@ SUBROUTINE dkes_printout(fz1p, fz1m, fz3p, fz3m, srces) c in the main DKES routine, only the STELLOPT routine changes this c and allocates these arrays IF (ALLOCATED(DKES_L11p) .and. (DKES_rad_dex > 0)) THEN -! PRINT *,'got here0' DKES_L11p(DKES_rad_dex) = L11p DKES_L33p(DKES_rad_dex) = L33p DKES_L31p(DKES_rad_dex) = L31p @@ -155,7 +154,6 @@ SUBROUTINE dkes_printout(fz1p, fz1m, fz3p, fz3m, srces) DKES_scal11(DKES_rad_dex) = scal11 DKES_scal33(DKES_rad_dex) = scal33 DKES_scal31(DKES_rad_dex) = scal13 -! PRINT *,'got here1' END IF c output results summary From 09baaa6524b5c25709b4bd67bc7328240449dd0a Mon Sep 17 00:00:00 2001 From: ajchcoelho Date: Wed, 13 Nov 2024 15:33:50 +0100 Subject: [PATCH 41/47] THRIFT: fix integration of DKES+PENTA in THRIFT --- PENTA/Sources/penta_interface_mod.f90 | 18 ++- THRIFT/Sources/thrift_boozer.f90 | 11 +- THRIFT/Sources/thrift_dkes.f90 | 40 +++-- THRIFT/Sources/thrift_funcs.f90 | 5 + THRIFT/Sources/thrift_penta.f90 | 221 ++++++++++++++++++++------ 5 files changed, 223 insertions(+), 72 deletions(-) diff --git a/PENTA/Sources/penta_interface_mod.f90 b/PENTA/Sources/penta_interface_mod.f90 index fd12a703c..ce43c7dbb 100644 --- a/PENTA/Sources/penta_interface_mod.f90 +++ b/PENTA/Sources/penta_interface_mod.f90 @@ -90,9 +90,12 @@ SUBROUTINE penta_set_ion_params(num_ion_in, Z_ion_in, miomp_in) num_ion_species = num_ion_in Z_ion_init = Z_ion_in miomp_init = miomp_in + ! Electrons + ions (added here because the function where this used to be defined [read_ion_params_namelist] is no longer called) AC 11/2024 + num_species = num_ion_species + 1_iknd RETURN END SUBROUTINE penta_set_ion_params + !! This subroutine is not called since in THRIFT we read the ions data from the h5 file SUBROUTINE read_penta_ion_params_namelist(filename,istat) USE safe_open_mod IMPLICIT NONE @@ -378,6 +381,17 @@ SUBROUTINE penta_deallocate_species IF (ALLOCATED(Er_test_vals)) DEALLOCATE(Er_test_vals) IF (ALLOCATED(QoT_i_vs_Er)) DEALLOCATE(QoT_i_vs_Er) IF (ALLOCATED(QoT_e_vs_Er)) DEALLOCATE(QoT_e_vs_Er) + IF (ALLOCATED(Gammas_ambi)) DEALLOCATE(Gammas_ambi) + IF (ALLOCATED(QoTs_ambi)) DEALLOCATE(QoTs_ambi) + IF (ALLOCATED(Flows_ambi)) DEALLOCATE(Flows_ambi) + IF (ALLOCATED(Jprl_ambi)) DEALLOCATE(Jprl_ambi) + IF (ALLOCATED(Jprl_parts)) DEALLOCATE(Jprl_parts) + IF (ALLOCATED(J_BS_ambi)) DEALLOCATE(J_BS_ambi) + IF (ALLOCATED(sigma_par_ambi)) DEALLOCATE(sigma_par_ambi) + IF (ALLOCATED(sigma_par_Spitzer_ambi)) DEALLOCATE(sigma_par_Spitzer_ambi) + IF (ALLOCATED(utor)) DEALLOCATE(utor) + IF (ALLOCATED(upol)) DEALLOCATE(upol) + RETURN END SUBROUTINE penta_deallocate_species @@ -520,8 +534,8 @@ SUBROUTINE penta_set_pprof(ne_in,dnedr_in,te_in,dtedr_in,ni_in,dnidr_in,ti_in,dt dtedr = dtedr_in ni = ni_in dnidr = dnidr_in - ti = ni_in - dtidr = dnidr_in + ti = ti_in + dtidr = dtidr_in RETURN END SUBROUTINE penta_set_pprof diff --git a/THRIFT/Sources/thrift_boozer.f90 b/THRIFT/Sources/thrift_boozer.f90 index 4921c7e3d..2688f6564 100644 --- a/THRIFT/Sources/thrift_boozer.f90 +++ b/THRIFT/Sources/thrift_boozer.f90 @@ -62,6 +62,10 @@ SUBROUTINE thrift_boozer(lscreen,iflag) CALL MPI_BCAST(mboz_xboozer,1,MPI_INTEGER,master,MPI_COMM_MYWORLD,ierr_mpi) CALL MPI_BCAST(nboz_xboozer,1,MPI_INTEGER,master,MPI_COMM_MYWORLD,ierr_mpi) CALL MPI_BCAST(ns_vmec,1,MPI_INTEGER,master,MPI_COMM_MYWORLD,ierr_mpi) + + !New -- introduced in order to have lsurf_boz calculated by master outside this routine (AC) + IF (myworkid /= master) ALLOCATE(lsurf_boz(ns_vmec)) + CALL MPI_BCAST(lsurf_boz,ns_vmec,MPI_LOGICAL,master,MPI_COMM_MYWORLD,ierr_mpi) #endif ! We need to setup the lsurf_boz array (this in now calculated outside this routine) !IF (ALLOCATED(lsurf_boz)) DEALLOCATE(lsurf_boz) @@ -229,8 +233,11 @@ SUBROUTINE thrift_boozer(lscreen,iflag) IF (ALLOCATED(beta_b)) DEALLOCATE(beta_b); ALLOCATE(beta_b(ns_b)); beta_b=0 IF (ALLOCATED(buco_b)) DEALLOCATE(buco_b); ALLOCATE(buco_b(ns_b)); buco_b=0 IF (ALLOCATED(bvco_b)) DEALLOCATE(bvco_b); ALLOCATE(bvco_b(ns_b)); bvco_b=0 - idx_b = 0 - idx_b(2:ns_b) = 1 + ! idx_b = 0 + ! idx_b(2:ns_b) = 1 + !NEW (comment lines above and introduce MERGE below) AC, 11/2024 + idx_b = MERGE(1,0,lsurf_boz) + ! iota_b = hiota pres_b = pres phip_b = phip diff --git a/THRIFT/Sources/thrift_dkes.f90 b/THRIFT/Sources/thrift_dkes.f90 index e3a38ebda..021820eec 100644 --- a/THRIFT/Sources/thrift_dkes.f90 +++ b/THRIFT/Sources/thrift_dkes.f90 @@ -45,8 +45,8 @@ SUBROUTINE thrift_dkes(lscreen,iflag) INTEGER, DIMENSION(:), ALLOCATABLE :: ik_dkes REAL(rprec), DIMENSION(:), ALLOCATABLE :: Earr_dkes, nuarr_dkes REAL(rprec), DIMENSION(:), POINTER :: f0p1, f0p2, f0m1, f0m2 - REAL(rprec) :: tcpu0, tcpu1, tcpui, tcput, tcpu, tcpua, dkes_efield,& - phi_temp + REAL(rprec) :: tcpu0, tcpu1, tcpui, tcput, tcpu, tcpua, & + phi_temp, stime, etime CHARACTER :: tb*1 CHARACTER*50 :: arg1(6) INTEGER :: numargs, iodata, iout_opt, idata, iout @@ -120,21 +120,22 @@ SUBROUTINE thrift_dkes(lscreen,iflag) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Note there are multile file write statements which can ! probably be removed so we only write what we need for PENTA + CALL second0(stime) DO k = mystart,myend ! Setup input arg1(1) = TRIM(proc_string) WRITE(arg1(2),'(i3)') ik_dkes(k) WRITE(arg1(3),'(e20.10)') nuarr_dkes(k) - WRITE(arg1(4),'(e20.10)') dkes_efield + WRITE(arg1(4),'(e20.10)') Earr_dkes(k) !dkes_efield arg1(5) = 'F' IF (lscreen .and. lfirst_pass) arg1(5) = 'T' WRITE(temp_str,'(i3.3)') k - arg1(6) = '_s' // TRIM(temp_str) + arg1(6) = '_k' // TRIM(temp_str) ier_phi = 0 ! We don't read the boozmn or wout file we've done that already - CALL dkes_input_prepare(arg1,6,dkes_input_file,ier_phi) - output_file= 'dkesout.' // TRIM(proc_string) // '_s' // TRIM(temp_str) - opt_file= 'opt_dkes.' // TRIM(proc_string) // '_s' // TRIM(temp_str) !DAS 2/21/2000 !Probably won't need - summary_file = 'results.' // TRIM(proc_string) //'_s' // TRIM(temp_str) !record file addition + CALL dkes_input_prepare_old(arg1,6,dkes_input_file,ier_phi) + output_file= 'dkesout.' // TRIM(proc_string) // '_k' // TRIM(temp_str) + opt_file= 'opt_dkes.' // TRIM(proc_string) // '_k' // TRIM(temp_str) !DAS 2/21/2000 !Probably won't need + summary_file = 'results.' // TRIM(proc_string) //'_k' // TRIM(temp_str) !record file addition ! OPEN INPUT AND OUTPUT FILES FOR READING (AND WRITING OUTPUT) idata = 7 iout = 30 @@ -244,11 +245,11 @@ SUBROUTINE thrift_dkes(lscreen,iflag) g31m, g13m, crs1m, crs3m) ENDIF ! This is a trick to get the arrays corretly sorted - DKES_rad_dex = i + DKES_rad_dex = k !i ---> I think w/ 'i' is completely wrong!! IF (.not. lfirst_pass) lscreen_dkes = .FALSE. CALL dkes_printout (f0p1, f0m1, f0p2, f0m2, srces0) - DKES_rad_dex = ik_dkes(i) - ! End trick + DKES_rad_dex = ik_dkes(k) !ik_dkes(i) -- same here !!! + ! End trickMPI_O CALL second0 (tcpu1); tcpu = tcpu1 - tcpu0; tcpu0 = tcpu1; tcput = tcput + tcpu; tcpua = tcput/i WRITE (ioout, 1100) tcpu CLOSE(unit=itab_out) @@ -261,7 +262,9 @@ SUBROUTINE thrift_dkes(lscreen,iflag) CLOSE(unit=ioout) CLOSE(unit=ioout_opt) lfirst_pass = .FALSE. + CALL second0(etime) END DO + PRINT *, 'I took ', etime-stime,'s.' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!! Parallel Work block Done !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -289,7 +292,7 @@ SUBROUTINE thrift_dkes(lscreen,iflag) CALL MPI_REDUCE(DKES_scal31, DKES_scal31, nruns_dkes, MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_MYWORLD, ierr_mpi) END IF #endif - IF (lscreen) WRITE(6,'(a)') ' ----------------- DKES CALCULATION (DONE) ----------------' + !IF (lscreen) WRITE(6,'(a)') ' ----------------- DKES CALCULATION (DONE) ----------------' IF (myworkid .ne. master) CALL read_wout_deallocate ! Sort out the DKES runs for later use @@ -299,16 +302,11 @@ SUBROUTINE thrift_dkes(lscreen,iflag) CALL mpialloc(DKES_D11, DKES_NK, DKES_NC, DKES_NE, myid_sharmem, 0, MPI_COMM_MYWORLD, win_dkes_d11) CALL mpialloc(DKES_D31, DKES_NK, DKES_NC, DKES_NE, myid_sharmem, 0, MPI_COMM_MYWORLD, win_dkes_d31) CALL mpialloc(DKES_D33, DKES_NK, DKES_NC, DKES_NE, myid_sharmem, 0, MPI_COMM_MYWORLD, win_dkes_d33) + IF (myworkid .eq. master) THEN - DO l = 1, nruns_dkes - i = MOD(l-1,DKES_NK)+1 - j = MOD(l-1,DKES_NK*DKES_NE) - j = FLOOR(REAL(j) / REAL(DKES_NK))+1 - k = CEILING(REAL(l) / REAL(DKES_NK*DKES_NE)) - DKES_D11(i,j,k) = (DKES_L11p(l) + DKES_L11m(l))*0.5 - DKES_D31(i,j,k) = (DKES_L31p(l) + DKES_L31m(l))*0.5 - DKES_D33(i,j,k) = (DKES_L33p(l) + DKES_L33m(l))*0.5 - END DO + DKES_D11 = RESHAPE( (DKES_L11p + DKES_L11m)*0.5, shape=(/DKES_NK, DKES_NC, DKES_NE/), order=(/2,3,1/) ) + DKES_D31 = RESHAPE( (DKES_L31p + DKES_L31m)*0.5, shape=(/DKES_NK, DKES_NC, DKES_NE/), order=(/2,3,1/) ) + DKES_D33 = RESHAPE( (DKES_L33p + DKES_L33m)*0.5, shape=(/DKES_NK, DKES_NC, DKES_NE/), order=(/2,3,1/) ) END IF END IF IF (lscreen) WRITE(6,'(a)') ' ------------------- DKES NEOCLASSICAL CALCULATION DONE ---------------------' diff --git a/THRIFT/Sources/thrift_funcs.f90 b/THRIFT/Sources/thrift_funcs.f90 index 8dc9df92e..c29abbcb6 100644 --- a/THRIFT/Sources/thrift_funcs.f90 +++ b/THRIFT/Sources/thrift_funcs.f90 @@ -64,6 +64,11 @@ SUBROUTINE update_vars() CALL get_prof_etapara(MIN(rho,SQRT(THRIFT_S(nsj-1))),mytime,THRIFT_ETAPARA(i,mytimestep)) ELSEIF(etapar_type == 'read_from_file') THEN CALL get_prof_eta(rho,mytime,THRIFT_ETAPARA(i,mytimestep)) + ELSEIF(etapar_type == 'dkespenta') THEN + IF(bootstrap_type /= 'dkespenta') THEN + STOP 'NOT POSSIBLE YET TO COMPUTE etapar FROM dkespenta WITHOUT BOOTSTRAP ALSO FROM dkespenta' + END IF + ! etapar will be computed when bootstrap code is called after update_vars ENDIF CALL get_prof_p(rho, mytime, THRIFT_P(i,mytimestep)) diff --git a/THRIFT/Sources/thrift_penta.f90 b/THRIFT/Sources/thrift_penta.f90 index 968fff8d2..bd04b2595 100644 --- a/THRIFT/Sources/thrift_penta.f90 +++ b/THRIFT/Sources/thrift_penta.f90 @@ -13,6 +13,7 @@ SUBROUTINE thrift_penta(lscreen,iflag) USE thrift_profiles_mod USE thrift_equil USE thrift_funcs + USE phys_const, ONLY: p_mass USE penta_interface_mod USE mpi_params USE mpi_inc @@ -28,16 +29,19 @@ SUBROUTINE thrift_penta(lscreen,iflag) ! Local Variables !----------------------------------------------------------------------- INTEGER :: ns_dkes, k, ier, j, i, ncstar, nestar, mystart, myend, & - mysurf - REAL(rprec) :: s, rho, iota, phip, chip, btheta, bzeta, bsq, vp, & - te, ne, dtedrho, dnedrho - REAL(rprec), DIMENSION(num_ion_species) :: ni,ti, dtidrho, dnidrho + mysurf, root_max_Er + REAL(rprec) :: s, rho + REAL(rprec), DIMENSION(:), ALLOCATABLE :: rho_k, iota, phip, chip, btheta, bzeta, bsq, vp, & + te, ne, dtedrho, dnedrho, EparB, JBS_PENTA, etapar_PENTA, rho_temp, J_temp, eta_temp + REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: ni,ti, dtidrho, dnidrho REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: D11, D13, D33 + TYPE(EZspline1_r8) :: J_spl, eta_spl + INTEGER :: bcs0(2) !----------------------------------------------------------------------- ! BEGIN SUBROUTINE !----------------------------------------------------------------------- IF (iflag < 0) RETURN - IF (lscreen) WRITE(6,'(a)') ' -------------------- NEOCLASSICAL FLUX USING PENTA -------------------' + IF (lscreen) WRITE(6,'(a)') ' -------------------- NEOCLASSICAL BOOTSTRAP USING PENTA -------------------' IF (lvmec) THEN ierr_mpi = 0 ! PENTA is parallelized over radial surfaces in this routine. @@ -47,52 +51,100 @@ SUBROUTINE thrift_penta(lscreen,iflag) END DO ! Break up work CALL MPI_CALC_MYRANGE(MPI_COMM_MYWORLD,1,ns_dkes,mystart,myend) - DO k = mystart,myend - ! Calc some information - mysurf = DKES_K(k) - s = DBLE(mysurf)/DBLE(ns_eq) - rho = sqrt(s) - ier = 0 - CALL EZSpline_interp(iota_spl, rho, iota, ier) - ier = 0 - CALL EZSpline_interp(phip_spl, rho, phip, ier) - ! PENTA wants d/ds note that this won't work if rho=0 - phip = 0.5*phip/rho - chip = iota * phip - ier = 0 - CALL EZSpline_interp(bu_spl, rho, btheta, ier) - ier = 0 - CALL EZSpline_interp(bv_spl, rho, bzeta, ier) - ier = 0 - CALL EZSpline_interp(bsq_spl, rho, bsq, ier) - ier = 0 - CALL EZSpline_interp(vp_spl, rho, vp, ier) - ! Vp = dV/dPHI need dVds with VMEC normalization - vp = vp * phip /(pi2*pi2) - ! Profiles - CALL get_prof_te(rho, THRIFT_T(mytimestep), te) - CALL get_prof_ne(rho, THRIFT_T(mytimestep), ne) - CALL get_prof_teprime(rho, THRIFT_T(mytimestep), dtedrho) - CALL get_prof_neprime(rho, THRIFT_T(mytimestep), dnedrho) - DO j = 1, num_ion_species - CALL get_prof_ti(rho, THRIFT_T(mytimestep), j, ti(j)) - CALL get_prof_ni(rho, THRIFT_T(mytimestep), j, ni(j)) - CALL get_prof_tiprime(rho, THRIFT_T(mytimestep), j, dtidrho(j)) - CALL get_prof_niprime(rho, THRIFT_T(mytimestep), j, dnidrho(j)) + + ALLOCATE(rho_k(ns_dkes),iota(ns_dkes),phip(ns_dkes),chip(ns_dkes),btheta(ns_dkes),bzeta(ns_dkes),bsq(ns_dkes),vp(ns_dkes),EparB(ns_dkes)) + ALLOCATE(te(ns_dkes),ne(ns_dkes),dtedrho(ns_dkes),dnedrho(ns_dkes)) + ALLOCATE(ni(ns_dkes,nion_prof),ti(ns_dkes,nion_prof),dtidrho(ns_dkes,nion_prof),dnidrho(ns_dkes,nion_prof)) + ALLOCATE(JBS_PENTA(ns_dkes),etapar_PENTA(ns_dkes)) + + JBS_PENTA = 0.0; etapar_PENTA = 0.0 + + IF (myworkid == master) THEN + DO k = 1, ns_dkes + mysurf = DKES_K(k) + s = DBLE(mysurf-1) / DBLE(ns_eq-1) + rho = SQRT(s) + rho_k(k) = rho + + ier = 0 + CALL EZSpline_interp(iota_spl, rho, iota(k), ier) + ier = 0 + CALL EZSpline_interp(phip_spl, rho, phip(k), ier) + ! PENTA wants d/ds note that this won't work if rho=0 + phip(k) = 0.5*phip(k)/rho + chip(k) = iota(k) * phip(k) + + ier = 0 + CALL EZSpline_interp(bu_spl, rho, btheta(k), ier) + ier = 0 + CALL EZSpline_interp(bv_spl, rho, bzeta(k), ier) + ier = 0 + CALL EZSpline_interp(bsq_spl, rho, bsq(k), ier) + ier = 0 + + CALL EZSpline_interp(vp_spl, rho, vp(k), ier) + ! Vp = dV/dPHI need dVds with VMEC normalization + vp(k) = vp(k) * phip(k) /(pi2*pi2) + + ! Profiles + CALL get_prof_te(rho, THRIFT_T(mytimestep), te(k)) + CALL get_prof_ne(rho, THRIFT_T(mytimestep), ne(k)) + CALL get_prof_teprime(rho, THRIFT_T(mytimestep), dtedrho(k)) + CALL get_prof_neprime(rho, THRIFT_T(mytimestep), dnedrho(k)) + DO j = 1, nion_prof + CALL get_prof_ti(rho, THRIFT_T(mytimestep), j, ti(k,j)) + CALL get_prof_ni(rho, THRIFT_T(mytimestep), j, ni(k,j)) + CALL get_prof_tiprime(rho, THRIFT_T(mytimestep), j, dtidrho(k,j)) + CALL get_prof_niprime(rho, THRIFT_T(mytimestep), j, dnidrho(k,j)) + END DO + + !EparB + EparB(k) = 0.0_rprec !! temporary... cannot do THRIFT_EPARB(mysurf,mytimestep) since this is defined in thrift grid, not vmec... interpolation? END DO - ! DKES Data - ncstar = COUNT(DKES_NUSTAR < 1E10) - nestar = COUNT(DKES_ERSTAR < 1E10) + + ! print *, 'master_thrift_penta: dkes_d11', DKES_D11 + END IF +#if defined(MPI_OPT) + CALL MPI_BCAST(rho_k,ns_dkes,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(iota,ns_dkes,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(phip,ns_dkes,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(btheta,ns_dkes,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(bzeta,ns_dkes,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(bsq,ns_dkes,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(vp,ns_dkes,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(te,ns_dkes,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(ne,ns_dkes,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(dtedrho,ns_dkes,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(dnedrho,ns_dkes,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(ti,ns_dkes*nion_prof,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(ni,ns_dkes*nion_prof,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(dtidrho,ns_dkes*nion_prof,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(dnidrho,ns_dkes*nion_prof,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(EparB,ns_dkes,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + ! VMEC quantities + CALL MPI_BCAST(eq_Aminor,1,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_BCAST(eq_Rmajor,1,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi) + +#endif + + ! DKES Data + ncstar = COUNT(DKES_NUSTAR < 1E10) + nestar = COUNT(DKES_ERSTAR < 1E10) + + DO k = mystart,myend ! PENTA - CALL PENTA_SET_ION_PARAMS(nion_prof, DBLE(Zatom_prof), Matom_prof/AMU) - CALL PENTA_SET_COMMANDLINE(-250.0_rprec,250.0_rprec,DKES_K(k),1,0.0_rprec,1,'','','') + CALL PENTA_SET_ION_PARAMS(nion_prof, DBLE(Zatom_prof), Matom_prof/p_mass) + + !!! WE NEED TO CHANGE THE EPARB --- it needs to interact with THRIFT...!!! + !!! ALSO: IS THIS -250:250 OK?? I PROPOSE TO MOVE THIS TO THE INPUT... + CALL PENTA_SET_COMMANDLINE(-250.0_rprec,250.0_rprec,DKES_K(k),1,EparB(k),1,'','','') CALL PENTA_ALLOCATE_SPECIES - CALL PENTA_SET_EQ_DATA(rho,eq_Aminor,eq_Rmajor,vp,chip,phip,iota,btheta,bzeta,bsq) - CALL PENTA_SET_PPROF(ne,dnedrho,te,dtedrho,ni,dnidrho,ti,dtidrho) - ! Check here we have D31 from DKES but need D13 in PENTA + CALL PENTA_SET_EQ_DATA(rho_k(k),eq_Aminor,eq_Rmajor,vp(k),chip(k),phip(k),iota(k),btheta(k),bzeta(k),bsq(k)) + CALL PENTA_SET_PPROF(ne(k),dnedrho(k)/eq_Aminor,te(k),dtedrho(k)/eq_Aminor,ni(k,:),dnidrho(k,:)/eq_Aminor,ti(k,:),dtidrho(k,:)/eq_Aminor) + ! MAKE CORRECTIONS ON D31 AND D33 -- values coming from DKES2 miss Bsq factors (see J. Lore documentation) CALL PENTA_SET_DKES_STAR(ncstar,nestar,DKES_NUSTAR(1:ncstar),DKES_ERSTAR(1:nestar), & - DKES_D11(mysurf,:,:),DKES_D31(mysurf,:,:),DKES_D33(mysurf,:,:)) + DKES_D11(k,:,:), DKES_D31(k,:,:)*SQRT(bsq(k)), DKES_D33(k,:,:)*bsq(k)) CALL PENTA_SET_BEAM(0.0_rprec) ! Zero becasue we don't read CALL PENTA_SET_U2() ! Leave blank for default value CALL PENTA_READ_INPUT_FILES(.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.) @@ -104,10 +156,85 @@ SUBROUTINE thrift_penta(lscreen,iflag) ! Now the basic steps CALL PENTA_RUN_2_EFIELD CALL PENTA_RUN_3_AMBIPOLAR + + ! Save JBS corresponding to the root that has the largest Er + ! This because whenever there are 2 stable roots, a rule of thumb is to pick the one with largest Er + root_max_Er = MAXLOC(Er_roots(1:num_roots),1) + JBS_PENTA(k) = J_BS_ambi(root_max_Er) + etapar_PENTA(k) = sigma_par_ambi(root_max_Er) + CALL PENTA_RUN_4_CLEANUP END DO + + !! Bootstrap interpolation onto THRIFT grid +#if defined(MPI_OPT) + CALL MPI_BARRIER(MPI_COMM_MYWORLD,ierr_mpi) + IF (myworkid == master) THEN + CALL MPI_REDUCE(MPI_IN_PLACE,JBS_PENTA,ns_dkes,MPI_DOUBLE_PRECISION,MPI_SUM,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_REDUCE(MPI_IN_PLACE,etapar_PENTA,ns_dkes,MPI_DOUBLE_PRECISION,MPI_SUM,master,MPI_COMM_MYWORLD,ierr_mpi) + ELSE + CALL MPI_REDUCE(JBS_PENTA,JBS_PENTA,ns_dkes,MPI_DOUBLE_PRECISION,MPI_SUM,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL MPI_REDUCE(etapar_PENTA,etapar_PENTA,ns_dkes,MPI_DOUBLE_PRECISION,MPI_SUM,master,MPI_COMM_MYWORLD,ierr_mpi) + DEALLOCATE(rho_k,iota,phip,chip,btheta,bzeta,bsq,vp,EparB) + DEALLOCATE(te,ne,dtedrho,dnedrho) + DEALLOCATE(ni,ti,dtidrho,dnidrho) + DEALLOCATE(JBS_PENTA,etapar_PENTA) + RETURN + ENDIF +#endif + + IF (myworkid == master) THEN + + ! Interpolate JBS_PENTA and etapar_PENTA at rho=0 and rho=1 + ALLOCATE(J_temp(ns_dkes+2),eta_temp(ns_dkes+2),rho_temp(ns_dkes+2)) + rho_temp(1) = 0.0 + rho_temp(2:ns_dkes+1) = rho_k + rho_temp(ns_dkes+2) = 1.0 + ! + J_temp(2:ns_dkes+1) = JBS_PENTA + J_temp(1) = J_temp(2) - (J_temp(3)-J_temp(2)) * rho_temp(2) / (rho_temp(3)-rho_temp(2)) + J_temp(ns_dkes+2) = JBS_PENTA(ns_dkes-1) + (JBS_PENTA(ns_dkes)-JBS_PENTA(ns_dkes-1)) * (1-rho_k(ns_dkes-1)) / (rho_k(ns_dkes)-rho_k(ns_dkes-1)) + ! + eta_temp(2:ns_dkes+1) = etapar_PENTA + eta_temp(1) = eta_temp(2) - (eta_temp(3)-eta_temp(2)) * rho_temp(2) / (rho_temp(3)-rho_temp(2)) + eta_temp(ns_dkes+2) = etapar_PENTA(ns_dkes-1) + (etapar_PENTA(ns_dkes)-etapar_PENTA(ns_dkes-1)) * (1-rho_k(ns_dkes-1)) / (rho_k(ns_dkes)-rho_k(ns_dkes-1)) + + ! Splines + bcs0=(/ 0, 0/) + !JBS + CALL EZspline_init(J_spl,ns_dkes+2,bcs0,ier) + J_spl%x1 = rho_temp + J_spl%isHermite = 1 + CALL EZspline_setup(J_spl,J_temp,ier,EXACT_DIM=.true.) + !etapar + CALL EZspline_init(eta_spl,ns_dkes+2,bcs0,ier) + eta_spl%x1 = rho_temp + eta_spl%isHermite = 1 + CALL EZspline_setup(eta_spl,eta_temp,ier,EXACT_DIM=.true.) + ! + DEALLOCATE(J_temp,eta_temp,rho_temp) + + ! Calculate J_BS and etapara in THRFIT GRID + DO i = 1, nsj + rho = SQRT( THRIFT_S(i) ) + CALL EZspline_interp(J_spl,rho,THRIFT_JBOOT(i,mytimestep),ier) + IF( etapar_type == 'dkespenta') CALL EZspline_interp(eta_spl,rho,THRIFT_ETAPARA(i,mytimestep),ier) + END DO + CALL EZspline_free(J_spl,ier) + CALL EZspline_free(eta_spl,ier) + + ! PRINT *, 'THRIFT_JBOOT=', THRIFT_JBOOT(:,mytimestep) + ! PRINT *, 'THRIFT_ETAPARA=', THRIFT_ETAPARA(:,mytimestep) + + DEALLOCATE(rho_k,iota,phip,chip,btheta,bzeta,bsq,vp,EparB) + DEALLOCATE(te,ne,dtedrho,dnedrho) + DEALLOCATE(ni,ti,dtidrho,dnidrho) + DEALLOCATE(JBS_PENTA,etapar_PENTA) + + END IF + ENDIF - IF (lscreen) WRITE(6,'(a)') ' ------------------- NEOCLASSICAL FLUX CALCULATION DONE ---------------------' + IF (lscreen) WRITE(6,'(a)') ' ------------------- NEOCLASSICAL BOOTSTRAP CALCULATION DONE ---------------------' RETURN !----------------------------------------------------------------------- ! END SUBROUTINE From a28e4efef190b3766b286cc62d2a14e8858c4f80 Mon Sep 17 00:00:00 2001 From: ajchcoelho Date: Wed, 13 Nov 2024 16:45:28 +0100 Subject: [PATCH 42/47] THRIFT: print to screen --- PENTA/Sources/penta_interface_mod.f90 | 9 ++++++--- THRIFT/Sources/thrift_dkes.f90 | 3 ++- THRIFT/Sources/thrift_penta.f90 | 2 ++ 3 files changed, 10 insertions(+), 4 deletions(-) diff --git a/PENTA/Sources/penta_interface_mod.f90 b/PENTA/Sources/penta_interface_mod.f90 index ce43c7dbb..e38c11c3a 100644 --- a/PENTA/Sources/penta_interface_mod.f90 +++ b/PENTA/Sources/penta_interface_mod.f90 @@ -1149,9 +1149,12 @@ SUBROUTINE penta_run_4_cleanup ! QQ write file with number of roots per surface! ! Write screen output - write(str_num,*) num_roots - write(*,'(f7.3,' // trim(adjustl(str_num)) // '(" ",e15.4))') & - roa_surf,er_roots(1:num_roots)/100._rknd + IF (lscreen) THEN + write(str_num,*) num_roots + write(*,'(f7.3,' // trim(adjustl(str_num)) // '(" ",e15.4))') & + roa_surf,er_roots(1:num_roots)/100._rknd + END IF + ! DEALLOCATE CALL penta_deallocate_species CALL penta_deallocate_dkescoeff diff --git a/THRIFT/Sources/thrift_dkes.f90 b/THRIFT/Sources/thrift_dkes.f90 index 021820eec..f3fd0edae 100644 --- a/THRIFT/Sources/thrift_dkes.f90 +++ b/THRIFT/Sources/thrift_dkes.f90 @@ -262,8 +262,9 @@ SUBROUTINE thrift_dkes(lscreen,iflag) CLOSE(unit=ioout) CLOSE(unit=ioout_opt) lfirst_pass = .FALSE. - CALL second0(etime) + IF(lscreen) WRITE(*, '(I0, A, I0, A, I0, A)') k, ' in [', mystart, ',', myend, '] completed' END DO + CALL second0(etime) PRINT *, 'I took ', etime-stime,'s.' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!! Parallel Work block Done diff --git a/THRIFT/Sources/thrift_penta.f90 b/THRIFT/Sources/thrift_penta.f90 index bd04b2595..a7d647b38 100644 --- a/THRIFT/Sources/thrift_penta.f90 +++ b/THRIFT/Sources/thrift_penta.f90 @@ -42,6 +42,8 @@ SUBROUTINE thrift_penta(lscreen,iflag) !----------------------------------------------------------------------- IF (iflag < 0) RETURN IF (lscreen) WRITE(6,'(a)') ' -------------------- NEOCLASSICAL BOOTSTRAP USING PENTA -------------------' + IF (lscreen) WRITE(6,'(A)') " /"," Er root(s) (V/cm)" + IF (lvmec) THEN ierr_mpi = 0 ! PENTA is parallelized over radial surfaces in this routine. From d7e6f6cb54d365d023acbb444675ac2cf5f6690a Mon Sep 17 00:00:00 2001 From: ajchcoelho Date: Wed, 13 Nov 2024 16:47:38 +0100 Subject: [PATCH 43/47] pySTEL: fix bug in plasma class --- pySTEL/libstell/plasma.py | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/pySTEL/libstell/plasma.py b/pySTEL/libstell/plasma.py index 1cab71de0..08380df7d 100644 --- a/pySTEL/libstell/plasma.py +++ b/pySTEL/libstell/plasma.py @@ -457,15 +457,15 @@ def write_THRIFT_profiles(self,time_array,filename=None): Te = np.tile(Te, (nt,1)).T # ni - ni = np.array( [self.get_density(ion,rho) for ion in self.ion_species] ) - ni = np.tile(ni[:,:,np.newaxis], (1,1,nt)).reshape((nrho,nt,self.num_ion_species)) + ni = np.repeat(ni[:,:,np.newaxis], nt, axis=2) + ni = np.transpose(ni, axes=[1,2,0]) # Ti - Ti = np.array( [self.get_temperature(ion,rho) for ion in self.ion_species] ) - Ti = np.tile(Ti[:,:,np.newaxis], (1,1,nt)).reshape((nrho,nt,self.num_ion_species)) - + Ti = np.repeat(Ti[:,:,np.newaxis], nt, axis=2) + Ti = np.transpose(Ti, axes=[1,2,0]) + hf = h5py.File(filename, 'w') hf.create_dataset('nrho', data=nrho) @@ -663,10 +663,11 @@ def print_SFINCS_list_namelist(self,roa_list,folder_path): file_content = f"""! Input file for SFINCS version 3. &general +RHSmode = 1 ambipolarSolve = .true. !.false. -Er_min = -10 -Er_max = 10 -NEr_ambipolarSolve = 10 +!Er_min = -10 +!Er_max = 10 +!NEr_ambipolarSolve = 10 / &geometryParameters @@ -677,7 +678,7 @@ def print_SFINCS_list_namelist(self,roa_list,folder_path): inputRadialCoordinateForGradients = 3 !the radial coordinate of the gradients given in species parameters is rN=sqrt(PHI/PHIEDGE) VMECRadialOption = 1 !get the nearest available flux surface from VMEC HALF grid -equilibriumFile = "wout_GIGA_v120.nc" +equilibriumFile = "wout_beta_2.nc" min_Bmn_to_load = 1e-4 / @@ -710,9 +711,9 @@ def print_SFINCS_list_namelist(self,roa_list,folder_path): &resolutionParameters Ntheta = 23 ! needs to be an odd number -Nzeta = 121 ! needs to be an odd number and at low collisionality might be needeed to be of the order 100 to converge +Nzeta = 91 ! needs to be an odd number and at low collisionality might be needeed to be of the order 100 to converge -Nxi = 100 +Nxi = 70 Nx = 6 solverTolerance = 1d-6 / From b6a4702b4bb0832bb74d67abc420e609b51c2b09 Mon Sep 17 00:00:00 2001 From: ajchcoelho Date: Wed, 13 Nov 2024 17:01:41 +0100 Subject: [PATCH 44/47] THRIFT: Er_min and Er_max now read from input namelist --- PENTA/Sources/penta_interface_mod.f90 | 6 ++++-- THRIFT/Sources/thrift_penta.f90 | 5 +---- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/PENTA/Sources/penta_interface_mod.f90 b/PENTA/Sources/penta_interface_mod.f90 index e38c11c3a..9fa249459 100644 --- a/PENTA/Sources/penta_interface_mod.f90 +++ b/PENTA/Sources/penta_interface_mod.f90 @@ -26,7 +26,7 @@ MODULE PENTA_INTERFACE_MOD iocheck, ie, ind_X, ind_A, ispec1, min_ind, iroot, num_roots REAL(rknd) :: Kmin, Kmax, epsabs, epsrel, Er_min, Er_max, B_Eprl, & U2, vth_e, loglambda, Er_test, abs_Er, min_Er, eaEr_o_kTe, cmin, & - cmax, emin, emax, sigma_par, sigma_par_Spitzer, J_BS + cmax, emin, emax, sigma_par, sigma_par_Spitzer, J_BS, Er_min_Vcm, Er_max_Vcm REAL(rknd), DIMENSION(NUM_ION_MAX) :: Z_ion_init, miomp_init REAL(rknd), DIMENSION(NUM_ROOTS_MAX) :: Er_roots REAL(rknd), DIMENSION(:), ALLOCATABLE :: ion_mass, Z_ion, vth_i, & @@ -49,7 +49,7 @@ MODULE PENTA_INTERFACE_MOD NAMELIST /run_params/ input_is_Er, log_interp, use_quanc8, & read_U2_file, Add_Spitzer_to_D33, num_Er_test, numKsteps, & kord_pprof, keord, kcord, Kmin, Kmax, epsabs, epsrel, Method, & - flux_cap, output_QoT_vs_Er, use_beam + flux_cap, output_QoT_vs_Er, use_beam, Er_min_Vcm, Er_max_Vcm !----------------------------------------------------------------------- ! SUBROUTINES @@ -79,6 +79,8 @@ SUBROUTINE init_penta_input sigma_par_Spitzer = 0.0E+0_rknd J_BS = 0.0E+0_rknd method = 'DKES' + Er_min_Vcm = -250.0_rknd + Er_max_Vcm = 250.0_rknd RETURN END SUBROUTINE init_penta_input diff --git a/THRIFT/Sources/thrift_penta.f90 b/THRIFT/Sources/thrift_penta.f90 index a7d647b38..3c1cb4a97 100644 --- a/THRIFT/Sources/thrift_penta.f90 +++ b/THRIFT/Sources/thrift_penta.f90 @@ -137,10 +137,7 @@ SUBROUTINE thrift_penta(lscreen,iflag) DO k = mystart,myend ! PENTA CALL PENTA_SET_ION_PARAMS(nion_prof, DBLE(Zatom_prof), Matom_prof/p_mass) - - !!! WE NEED TO CHANGE THE EPARB --- it needs to interact with THRIFT...!!! - !!! ALSO: IS THIS -250:250 OK?? I PROPOSE TO MOVE THIS TO THE INPUT... - CALL PENTA_SET_COMMANDLINE(-250.0_rprec,250.0_rprec,DKES_K(k),1,EparB(k),1,'','','') + CALL PENTA_SET_COMMANDLINE(Er_min_Vcm,Er_max_Vcm,DKES_K(k),1,EparB(k),1,'','','') CALL PENTA_ALLOCATE_SPECIES CALL PENTA_SET_EQ_DATA(rho_k(k),eq_Aminor,eq_Rmajor,vp(k),chip(k),phip(k),iota(k),btheta(k),bzeta(k),bsq(k)) CALL PENTA_SET_PPROF(ne(k),dnedrho(k)/eq_Aminor,te(k),dtedrho(k)/eq_Aminor,ni(k,:),dnidrho(k,:)/eq_Aminor,ti(k,:),dtidrho(k,:)/eq_Aminor) From cef7043d1a37174c612cc5da05a40383c097993b Mon Sep 17 00:00:00 2001 From: ajchcoelho Date: Wed, 13 Nov 2024 17:02:48 +0100 Subject: [PATCH 45/47] THRIFT: Close PENTA files in interface --- PENTA/Sources/penta_interface_mod.f90 | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/PENTA/Sources/penta_interface_mod.f90 b/PENTA/Sources/penta_interface_mod.f90 index 9fa249459..aafbea630 100644 --- a/PENTA/Sources/penta_interface_mod.f90 +++ b/PENTA/Sources/penta_interface_mod.f90 @@ -1160,6 +1160,17 @@ SUBROUTINE penta_run_4_cleanup ! DEALLOCATE CALL penta_deallocate_species CALL penta_deallocate_dkescoeff + + ! Close files + ! Close output files + Close(iu_flux_out) + Close(iu_pprof_out) + Close(iu_fvEr_out) + Close(iu_QoTvEr_out) + Close(iu_flows_out) + Close(iu_flowvEr_out) + Close(iu_Jprl_out) + Close(iu_contraflows_out) END SUBROUTINE penta_run_4_cleanup From f65ce101947877ad0590a9c23bde8866d880652a Mon Sep 17 00:00:00 2001 From: ajchcoelho Date: Wed, 13 Nov 2024 17:17:08 +0100 Subject: [PATCH 46/47] PENTA: comment lscreen --- PENTA/Sources/penta_interface_mod.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PENTA/Sources/penta_interface_mod.f90 b/PENTA/Sources/penta_interface_mod.f90 index aafbea630..de4d88a0f 100644 --- a/PENTA/Sources/penta_interface_mod.f90 +++ b/PENTA/Sources/penta_interface_mod.f90 @@ -1151,11 +1151,11 @@ SUBROUTINE penta_run_4_cleanup ! QQ write file with number of roots per surface! ! Write screen output - IF (lscreen) THEN + !IF (lscreen) THEN write(str_num,*) num_roots write(*,'(f7.3,' // trim(adjustl(str_num)) // '(" ",e15.4))') & roa_surf,er_roots(1:num_roots)/100._rknd - END IF + !END IF ! DEALLOCATE CALL penta_deallocate_species From e81300ae5b388f74b1604170e434de8326819815 Mon Sep 17 00:00:00 2001 From: ajchcoelho Date: Thu, 14 Nov 2024 11:10:37 +0100 Subject: [PATCH 47/47] THRIFT: Improve prints to command line --- PENTA/Sources/penta_interface_mod.f90 | 4 ++-- THRIFT/Sources/thrift_dkes.f90 | 6 +++--- THRIFT/Sources/thrift_penta.f90 | 3 ++- 3 files changed, 7 insertions(+), 6 deletions(-) diff --git a/PENTA/Sources/penta_interface_mod.f90 b/PENTA/Sources/penta_interface_mod.f90 index de4d88a0f..27182099d 100644 --- a/PENTA/Sources/penta_interface_mod.f90 +++ b/PENTA/Sources/penta_interface_mod.f90 @@ -861,8 +861,8 @@ SUBROUTINE penta_run_2_efield Else Er_test_vals(min_ind) = Er_test_vals(min_ind + 1)/2._rknd EndIf - Write(*,'(a,i4,a,f10.3)') 'Cannot use Er=0 with log_interp, using Er(', & - min_ind, ') = ', Er_test_vals(min_ind) + ! Write(*,'(a,i4,a,f10.3)') 'Cannot use Er=0 with log_interp, using Er(', & + ! min_ind, ') = ', Er_test_vals(min_ind) EndIf ! Loop over Er to get fluxes as a function of Er Do ie = 1,num_Er_test diff --git a/THRIFT/Sources/thrift_dkes.f90 b/THRIFT/Sources/thrift_dkes.f90 index f3fd0edae..a0c91b901 100644 --- a/THRIFT/Sources/thrift_dkes.f90 +++ b/THRIFT/Sources/thrift_dkes.f90 @@ -124,12 +124,12 @@ SUBROUTINE thrift_dkes(lscreen,iflag) DO k = mystart,myend ! Setup input arg1(1) = TRIM(proc_string) - WRITE(arg1(2),'(i3)') ik_dkes(k) + WRITE(arg1(2),'(i4)') ik_dkes(k) WRITE(arg1(3),'(e20.10)') nuarr_dkes(k) WRITE(arg1(4),'(e20.10)') Earr_dkes(k) !dkes_efield arg1(5) = 'F' IF (lscreen .and. lfirst_pass) arg1(5) = 'T' - WRITE(temp_str,'(i3.3)') k + WRITE(temp_str,'(i4.4)') k arg1(6) = '_k' // TRIM(temp_str) ier_phi = 0 ! We don't read the boozmn or wout file we've done that already CALL dkes_input_prepare_old(arg1,6,dkes_input_file,ier_phi) @@ -265,7 +265,7 @@ SUBROUTINE thrift_dkes(lscreen,iflag) IF(lscreen) WRITE(*, '(I0, A, I0, A, I0, A)') k, ' in [', mystart, ',', myend, '] completed' END DO CALL second0(etime) - PRINT *, 'I took ', etime-stime,'s.' + WRITE(*, '(A, I0, A, F8.2, A, I0, A)') 'Processor ', myworkid, ' took ', etime-stime,'s to compute ', myend-mystart+1, ' DKES coefficients.' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!! Parallel Work block Done !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/THRIFT/Sources/thrift_penta.f90 b/THRIFT/Sources/thrift_penta.f90 index 3c1cb4a97..4eeb18662 100644 --- a/THRIFT/Sources/thrift_penta.f90 +++ b/THRIFT/Sources/thrift_penta.f90 @@ -42,7 +42,7 @@ SUBROUTINE thrift_penta(lscreen,iflag) !----------------------------------------------------------------------- IF (iflag < 0) RETURN IF (lscreen) WRITE(6,'(a)') ' -------------------- NEOCLASSICAL BOOTSTRAP USING PENTA -------------------' - IF (lscreen) WRITE(6,'(A)') " /"," Er root(s) (V/cm)" + IF (lscreen) Write(*,*) " /"," Er root(s) (V/cm)" IF (lvmec) THEN ierr_mpi = 0 @@ -174,6 +174,7 @@ SUBROUTINE thrift_penta(lscreen,iflag) ELSE CALL MPI_REDUCE(JBS_PENTA,JBS_PENTA,ns_dkes,MPI_DOUBLE_PRECISION,MPI_SUM,master,MPI_COMM_MYWORLD,ierr_mpi) CALL MPI_REDUCE(etapar_PENTA,etapar_PENTA,ns_dkes,MPI_DOUBLE_PRECISION,MPI_SUM,master,MPI_COMM_MYWORLD,ierr_mpi) + CALL FLUSH(6) DEALLOCATE(rho_k,iota,phip,chip,btheta,bzeta,bsq,vp,EparB) DEALLOCATE(te,ne,dtedrho,dnedrho) DEALLOCATE(ni,ti,dtidrho,dnidrho)