From c5a4b180bc2a9280ab9e509ef7820193d66d53e5 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Mon, 29 Jan 2024 01:08:28 +0100 Subject: [PATCH] src/gray_params.f90: replace magic numbers with enums 1. Introduces enumerations (and some booleans) intended to replace all the magic numbers used throughout the code to represent multiple choices. 2. Replace the gray_params.sh script a new one that automatically generates code for all the GRAY parameters by parsing gray_params.f90. 3. Also generate extra code to accept the enum identifiers as valid values in the configuration files and command line arguments. 4. Set sensible default values for parameters that are rarely changes. --- Makefile | 6 +- input/gray.ini | 118 +++++++++--------- src/equilibrium.f90 | 10 +- src/gray_core.f90 | 9 +- src/gray_params.awk | 69 +++++++++++ src/gray_params.f90 | 196 ++++++++++++++++++++++-------- src/gray_params.sh | 30 ----- src/multipass.f90 | 5 +- src/pec.f90 | 2 +- tests/06-ITER-startup/__init__.py | 2 +- 10 files changed, 291 insertions(+), 156 deletions(-) create mode 100644 src/gray_params.awk delete mode 100644 src/gray_params.sh diff --git a/Makefile b/Makefile index 4eb8947..e49dd70 100644 --- a/Makefile +++ b/Makefile @@ -197,9 +197,9 @@ $(OBJDIR)/%.d: $(SRCDIR)/%.f90 | $(OBJDIR) $(OBJDIR)/%.o: $(SRCDIR)/%.f90 | $(OBJDIR) $(INCDIR) $(FC) -cpp $(CPPFLAGS) $(FFLAGS) -c '$<' -o '$@' -# Generate Fortran code from shell script -$(INCDIR)/%.inc: $(SRCDIR)/%.sh | $(INCDIR) - sh '$<' > '$@' +# Generate Fortran code from awk script +$(INCDIR)/%.inc: $(SRCDIR)/%.awk | $(INCDIR) + awk -f '$<' > '$@' # Create directories $(DIRS): diff --git a/input/gray.ini b/input/gray.ini index 60a0aee..8acec10 100644 --- a/input/gray.ini +++ b/input/gray.ini @@ -3,7 +3,7 @@ nrayr = 1 nrayth = 16 ; Normalized maximum radius of the beam power -rwmax = 1 +rwmax = 1.0 ; Switch between simple raytracing (0) and beamtracing (1) igrad = 0 @@ -13,29 +13,28 @@ igrad = 0 ; G-EQDSK file; when negative on a simple limiter at R=`rwall`, see below. ipass = 1 -; Whether to compute the wave polarisation -ipol = 0 +; Whether to compute the wave polarisation (from chi, psi) +ipol = .false. ; Step size (cm) for the numerical integration dst = 0.1 ; Max number of integration steps nstep = 12000 ; Choice of the integration variable -; 0: path length (s) -; 1: "time" (actually c⋅t) -; 2: real part of the eikonal (S_r) -idst = 0 +; STEP_ARCLEN: arc length (s) +; STEP_TIME: time (actually c⋅t) +; STEP_PHASE: phase (actually real part of eikonal S_r=k₀⋅φ) +idst = STEP_ARCLEN [ecrh_cd] ; Choice of the power absorption model -; 0: no absorption at all -; 1: weakly relativistic -; 2: fully relativistic (faster variant) -; 3: fully relativistic (slower variant) -; 4: tenuous plasma (very fast) -; Note: iwarm>0 is required for current driver -iwarm = 2 +; ABSORP_OFF: no absorption at all +; ABSORP_WEAK: weakly relativistic +; ABSORP_FULL: fully relativistic (faster variant) +; ABSORP_FULL_ALT: fully relativistic (slower variant) +; Note: iwarm /= ABSORP_OFF is required for current drive +iwarm = ABSORP_FULL ; Order of the electron Larmor radius expansion ; (used by some absorption models) @@ -48,11 +47,11 @@ ilarm = 5 imx = -20 ; Current drive model -; 0: no current drive at all -; 1: Cohen -; 2: no trapping -; 3: Neoclassical -ieccd = 3 +; CD_OFF: no current drive at all +; CD_COHEN: Cohen +; CD_NO_TRAP: no trapping +; CD_NEOCLASSIC: Neoclassical +ieccd = CD_NEOCLASSIC [antenna] @@ -61,22 +60,22 @@ alpha = 45 ; Poloidal angle (positive → up) beta = 0 ; Toroidal angle (positive → right) ; Injected power (MW) -power = 1 +power = 1.0 ; Polarisation mode -; 1: ordinary (O) -; 2: extraordinary (X) -iox = 1 +; MODE_O: ordinary (O) +; MODE_X: extraordinary (X) +iox = MODE_X ; Alternatively, parameters of the polarisation ellipse ; χ: angle between the principal axes and the (x,y) axes ; ψ: atan(ε), where ε is the ellipticity chi = 0 psi = 0 -; Beam description kind -; 0: simple beam shape -; 1: 1D table -; 2: 2D table +; Beam parameters format +; BEAM_0D: fixed beam parameters +; BEAM_1D: 1D steering angle table +; BEAM_2D: 2D steering angles table ibeam = 0 ; Filepath of the beam data (relative to this file) @@ -85,40 +84,40 @@ filenm = "beamdata.txt" [equilibrium] ; MHD equilibrium kind -; 0: vacuum (i.e. no plasma at all) -; 1: analytical -; 2: G-EQDSK format - data valid on the whole domain -; 3: G-EQDSK format - data valid only inside the LCFS -iequil = 3 +; EQ_VACUUM: vacuum (i.e. no plasma at all) +; EQ_ANALYTICAL: analytical model +; EQ_EQDSK_FULL: G-EQDSK format - data valid on the whole domain +; EQ_EQDSK_PARTIAL: G-EQDSK format - data valid only inside the LCFS +iequil = EQ_EQDSK_FULL ; Filepath of the equilibrium data (relative to this file) filenm = "magneticdata.eqdsk" ; COCOS index icocos = 0 -; Normalisation of the poloidal function -; 0: G-EQDSK, ψ → |ψ - ψ(edge)|/|ψ(axis) - ψ(edge)| -; 1: no normalisation -ipsinorm = 0 +; Whether the poloidal function is normalised (G-EQDSK) +; false: is not normalised, ψ → |ψ - ψ(edge)|/|ψ(axis) - ψ(edge)| +; true: is already normalised +ipsinorm = false ; G-EQDSK format parameters ; Whether header starts with a description, a.k.a identification string -idesc = 1 +idesc = true ; Whether the records have variable length ; Note: some non-compliant programs output numbers formatted with variable length ; instead of using the single (5e16.9) specifier. -ifreefmt = 0 +ifreefmt = false ; Position of the X point -; -1: bottom -; 0: no X point -; +1: top -ixp = 0 +; X_IS_MISSING: No X point +; X_AT_TOP: At the top of the plasma +; X_AT_BOTTOM: At the bottom of the plasma +ixp = X_IS_MISSING ; Tension of splines ; Note: 0 means perfect interpolation ssplps = 0.005 ; for ψ(R,Z), normalised poloidal flux -ssplf = 0.01 ; for I(ψ)=R⋅B_T, poloidal current function +ssplf = 0.01 ; for F(ψ)=R⋅B_T, poloidal current function ; Sign of toroidal field/current (used when COCOS index is 0,10) ; When viewing from above: +1 → counter-clockwise, -1 → clockwise @@ -131,16 +130,16 @@ factb = 1 [profiles] ; (input) plasma profiles parameters -; Profiles kind -; 0: analytical -; 1: numerical -iprof = 1 +; Plasma profiles kind +; PROF_ANALYTIC: analytical model +; PROF_NUMERIC: numerical data (ρ, n_e, T_e, table) +iprof = PROF_NUMERIC -; Profile radial coordinate -; 0: ρ_t = √Φ (where Φ is the normalised toroidal flux) -; 1: ρ_p = √ψ (where ψ is the normalised poloidal flux) -; 2: normalised poloidal flux ψ -irho = 0 +; Plasma profiles radial coordinate +; RHO_TOR: ρ_t = √Φ (where Φ is the normalised toroidal flux) +; RHO_POL: ρ_p = √ψ (where ψ is the normalised poloidal flux) +; RHO_PSI: normalised poloidal flux ψ +irho = RHO_TOR ; Filepath of the equilibrium (relative to this file) filenm = "profiles.txt" @@ -153,11 +152,12 @@ sspld = 0.1 factte = 1 factne = 1 -; Choice of model for rescaling the temperature/density with the -; magnetic field -; 1: costant Greenwald density (ab=1) -; 2: no rescaling (ab=0) -iscal = 2 +; Choice of model for rescaling the temperature/density +; with the magnetic field (if factb ≠ 0) +; SCALE_OFF: don't rescale at all +; SCALE_COLLISION: scale while preserving collisionality +; SCALE_GREENWALD: scale while preserving the Greenwald fraction +iscal = SCALE_OFF [output] @@ -165,8 +165,8 @@ iscal = 2 ; ECRH&CD profiles grid: ; Radial coordinate -; 1: ρ_p = √ψ (where ψ is the normalised poloidal flux) -; 2: ρ_t = √Φ (where Φ is the normalised toroidal flux) +; RHO_TOR: ρ_t = √Φ (where Φ is the normalised toroidal flux) +; RHO_POL: ρ_p = √ψ (where ψ is the normalised poloidal flux) ipec = 1 ; Number of points nrho = 501 diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index e45eebd..ee51ae8 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -114,7 +114,7 @@ contains end if ! get size of main arrays and allocate them - if (params%idesc == 1) then + if (params%idesc) then read (u,'(a48,3i4)') string,idum,nr,nz else read (u,*) nr, nz @@ -132,7 +132,7 @@ contains data%qpsi(nr)) ! Store 0D data and main arrays - if (params%ifreefmt==1) then + if (params%ifreefmt) then read (u, *) dr, dz, data%rvac, rleft, zmid read (u, *) data%rax, data%zax, psiaxis, psiedge, xdum read (u, *) xdum, xdum, xdum, xdum, xdum @@ -166,7 +166,7 @@ contains ! Load plasma boundary data if(nbnd > 0) then allocate(data%rbnd(nbnd), data%zbnd(nbnd)) - if (params%ifreefmt == 1) then + if (params%ifreefmt) then read(u, *) (data%rbnd(i), data%zbnd(i), i=1,nbnd) else read(u, '(5e16.9)') (data%rbnd(i), data%zbnd(i), i=1,nbnd) @@ -176,7 +176,7 @@ contains ! Load limiter data if(nlim > 0) then allocate(data%rlim(nlim), data%zlim(nlim)) - if (params%ifreefmt == 1) then + if (params%ifreefmt) then read(u, *) (data%rlim(i), data%zlim(i), i=1,nlim) else read(u, '(5e16.9)') (data%rlim(i), data%zlim(i), i=1,nlim) @@ -201,7 +201,7 @@ contains ! Normalize psin data%psia = psiedge - psiaxis - if(params%ipsinorm == 0) data%psin = (data%psin - psiaxis)/data%psia + if(.not. params%ipsinorm) data%psin = (data%psin - psiaxis)/data%psia end subroutine read_eqdsk diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 54ebb26..5e100b0 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -150,7 +150,7 @@ contains iroff,yynext,yypnext,yw0,ypw0, & stnext,p0ray,taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv) - if(params%raytracing%ipol == 0) then + if(.not. params%raytracing%ipol) then if(params%antenna%iox == 2) then ! only X mode on 1st pass cpl0 = [zero, one] else ! only O mode on 1st pass @@ -305,7 +305,7 @@ contains if(iop(jk) == 1 .and. ip==1) then ! * 1st entrance on 1st pass (ray hasn't entered in plasma yet) => continue current pass - if(params%raytracing%ipol == 0) then ! + IF single mode propagation + if(.not. params%raytracing%ipol) then ! + IF single mode propagation cpl = cpl0 p0ray(jk) = p0ray(jk)*cpl(iox) else if(cpl(iox) < etaucr) then ! + ELSE IF low coupled power for current mode => de-activate derived rays @@ -1842,7 +1842,8 @@ contains ! subroutine arguments real(wp_), dimension(6, nray), intent(in) :: ywrk0 real(wp_), intent(in) :: bres - integer, intent(in) :: sox, ipol + integer, intent(in) :: sox + logical, intent(in) :: ipol real(wp_), intent(inout) :: psipol0, chipol0 complex(wp_), dimension(nray), intent(out) :: ext0, eyt0 @@ -1860,7 +1861,7 @@ contains k=1 end if - if(ipol == 0) then + if(.not. ipol) then xmv=ywrk0(1:3,jk)*0.01_wp_ ! convert from cm to m anv=ywrk0(4:6,jk) rm=sqrt(xmv(1)**2+xmv(2)**2) diff --git a/src/gray_params.awk b/src/gray_params.awk new file mode 100644 index 0000000..026de74 --- /dev/null +++ b/src/gray_params.awk @@ -0,0 +1,69 @@ +#/usr/bin/env awk -f + +# This script generates the very repetitive code for parsing the +# GRAY parameters from a string. The output is embedded +# into gray_params.f90 using the #include cpp directive. + +BEGIN { + # set file to open + ARGV[1] = "src/gray_params.f90" + ARGC = 2 + + # start enum conversion switch + print "select case (value)" +} + +# match an enum constant +/enumerator :: [A-Z0-9_]+/ { + # generate conversion case + print " case ('"$3"')" + print " temp = '"$5"'" +} + +# match start of a parameters set +/type [[:alnum:]_]+_parameters$/ { + if (!done_enum) { + # close enum switch + print " case default" + print " temp = value" + print "end select\n" + # start params switch + print "select case (name)" + done_enum = 1 + } + # set the current set + split($2, words, "_parameters") + if (words[1] != "gray") set = words[1] + # change separator for the next step + FS = "::" +} + +# match a parameter +/ :: [[:alnum:]_]+/ { + if (!set) next + # strip comments + sub(/!.*$/, "") + for (i = 1; i <= split($2, ids, ","); i++) { + # extract identifier + match(ids[i], /^ *[[:alnum:]_]+/) + id = substr(ids[i], RSTART + 1, RLENGTH - 1) + # generate parameter case + print " case ('"set"."id"')" + print " read (temp, *, iostat=err) params%"set"%"id + print " if (err > 0) err = ERR_VALUE" + } +} + +# match end of a parameters set +/end type$/ { + set = "" + FS = " " +} + +END { + # close the parameters switch + print " case default" + print " ! unknown parameter" + print " err = ERR_UNKNOWN" + print "end select" +} diff --git a/src/gray_params.f90 b/src/gray_params.f90 index 537fc44..0cdf4c3 100644 --- a/src/gray_params.f90 +++ b/src/gray_params.f90 @@ -7,15 +7,96 @@ module gray_params integer, parameter :: lenfnm = 256 integer, parameter :: headw = 132, headl = 21 + ! Polarisation mode + enum, bind(c) + enumerator :: polar_enum = -1 + enumerator :: MODE_O = 1 ! ordinary mode (O) + enumerator :: MODE_X = 2 ! extraordinary mode (X) + end enum + + ! Beam parameters format + enum, bind(c) + enumerator :: beam_enum = -1 + enumerator :: BEAM_0D = 0 ! fixed beam parameters + enumerator :: BEAM_1D = 1 ! 1D steering angle table + enumerator :: BEAM_2D = 2 ! 2D steering angles table + end enum + + ! Position of the X point + enum, bind(c) + enumerator :: x_point_enum = -1 + enumerator :: X_IS_MISSING = 0 ! no X point + enumerator :: X_AT_TOP = +1 ! at the top of the plasma + enumerator :: X_AT_BOTTOM = -1 ! at the bottom of the plasma + end enum + + ! MHD equilibrium kind + enum, bind(c) + enumerator :: equil_enum = -1 + enumerator :: EQ_VACUUM = 0 ! vacuum (i.e. no plasma at all) + enumerator :: EQ_ANALYTICAL = 1 ! analytical model + enumerator :: EQ_EQDSK_FULL = 2 ! G-EQDSK format - data valid on the whole domain + enumerator :: EQ_EQDSK_PARTIAL = 3 ! G-EQDSK format - data valid only inside the LCFS + end enum + + ! Temperature/density scaling model + enum, bind(c) + enumerator :: scaling_enum = -1 + enumerator :: SCALE_COLLISION = 0 ! preserve collisionality + enumerator :: SCALE_GREENWALD = 1 ! preserve Greenwald fraction + enumerator :: SCALE_OFF = 2 ! don't rescale at all + end enum + + ! Radial profile coordinate + enum, bind(c) + enumerator :: rho_enum = -1 + enumerator :: RHO_TOR = 0 ! ρ_t = √Φ (where Φ is the normalised toroidal flux) + enumerator :: RHO_POL = 1 ! ρ_p = √ψ (where ψ is the normalised poloidal flux) + enumerator :: RHO_PSI = 2 ! normalised poloidal flux ψ + end enum + + ! Plasma profiles kind + enum, bind(c) + enumerator :: profiles_enum = -1 + enumerator :: PROF_ANALYTIC = 0 ! analytical model + enumerator :: PROF_NUMERIC = 1 ! numerical data (ρ, n_e, T_e, table) + end enum + + ! Choice of the integration variable + enum, bind(c) + enumerator :: step_enum = -1 + enumerator :: STEP_ARCLEN = 0 ! arc length (s) + enumerator :: STEP_TIME = 1 ! time (actually c⋅t) + enumerator :: STEP_PHASE = 2 ! phase (actually real part of eikonal S_r=k₀⋅φ) + end enum + + ! Absorption model + enum, bind(c) + enumerator :: absorption_enum = -1 + enumerator :: ABSORP_OFF = 0 ! no absorption at all + enumerator :: ABSORP_WEAK = 1 ! weakly relativistic + enumerator :: ABSORP_FULL = 2 ! fully relativistic (faster variant) + enumerator :: ABSORP_FULL_ALT = 3 ! fully relativistic (slower variant) + end enum + + ! Current drive model + enum, bind(c) + enumerator :: current_drive_enum = -1 + enumerator :: CD_OFF = 0 ! no current drive at all + enumerator :: CD_COHEN = 1 ! Cohen + enumerator :: CD_NO_TRAP = 2 ! no trapping + enumerator :: CD_NEOCLASSIC = 3 ! Neoclassical + end enum + ! Antenna/wave launcher parameters type antenna_parameters ! From gray_params.data: - real(wp_) :: alpha, beta ! Launching angles - real(wp_) :: power ! Initial power - real(wp_) :: psi, chi ! Initial polarisation ellipse parameters - integer :: iox ! Initial wave mode - integer :: ibeam ! Beam kind - character(len=lenfnm) :: filenm ! beamdata.txt filename + real(wp_) :: alpha, beta ! Launching angles + real(wp_) :: power = 1.0_wp_ ! Initial power + real(wp_) :: psi = 0, chi = 0 ! Initial polarisation ellipse parameters + integer(kind(polar_enum)) :: iox ! Initial wave mode + integer(kind(beam_enum)) :: ibeam ! Beam parameters format + character(len=lenfnm) :: filenm ! beamdata.txt filename ! From beamdata.txt: real(wp_) :: fghz ! EC frequency @@ -27,54 +108,57 @@ module gray_params ! MHD equilibrium parameters type equilibrium_parameters - real(wp_) :: ssplps, ssplf ! Spline tensions for ψ(R,Z), I(ψ) - real(wp_) :: factb ! Rescaling factor for B - integer :: sgnb, sgni ! Sign of B, I - integer :: ixp ! Position of X point - integer :: iequil ! Equilibrium kind - integer :: icocos ! COCOS index - integer :: ipsinorm ! Normalisation of ψ - integer :: idesc, ifreefmt ! G-EQDSK quirks - character(len=lenfnm) :: filenm ! Equilibrium data filepath + real(wp_) :: ssplps = 0.005_wp_ ! Spline tension for ψ(R,Z) + real(wp_) :: ssplf = 0.01_wp_ ! Spline tension for F(ψ) + real(wp_) :: factb = 1.0_wp_ ! Rescaling factor for B + integer :: sgnb = 0, sgni = 0 ! Sign of B, I + integer(kind(x_point_enum)) :: ixp = X_IS_MISSING ! Position of X point + integer(kind(equil_enum)) :: iequil = EQ_EQDSK_FULL ! Equilibrium kind + integer :: icocos = 3 ! COCOS index + logical :: ipsinorm = .false. ! Whether ψ is normalised + logical :: idesc = .true. ! G-EQDSK quirks + logical :: ifreefmt = .false. ! + character(len=lenfnm) :: filenm ! Equilibrium data filepath end type ! Kinetic plasma profiles type profiles_parameters - real(wp_) :: psnbnd ! Plasma boundary (ψ: ne(ψ)=0) - real(wp_) :: sspld ! Density spline tension - real(wp_) :: factne, factte ! Rescale factors for n_e, T_e - integer :: iscal ! Rescaling model - integer :: irho ! Radial coordinate - integer :: iprof ! Profile kind - character(len=lenfnm) :: filenm ! Profiles data filepath + real(wp_) :: psnbnd = 1.01_wp_ ! Plasma boundary (ψ: ne(ψ)=0) + real(wp_) :: sspld = 0.1_wp_ ! Density spline tension + real(wp_) :: factne = 1.0_wp_ ! Rescale factor for n_e + real(wp_) :: factte = 1.0_wp_ ! Rescale factor for T_e + integer(kind(scaling_enum)) :: iscal = SCALE_OFF ! Rescaling model for n_e,T_e + integer(kind(rho_enum)) :: irho ! Radial coordinate + integer(kind(profiles_enum)) :: iprof ! Profile kind + character(len=lenfnm) :: filenm ! Profiles data filepath end type ! Raytracing parameters type raytracing_parameters - real(wp_) :: dst ! Integration step size - real(wp_) :: rwmax ! Normalized maximum radius of beam power - integer :: nray, nrayr, nrayth ! Total/radial/angular number of rays - integer :: igrad ! Complex eikonal switch - integer :: nstep ! Max number of integration steps - integer :: idst ! Choice of the integration variable - integer :: ipass ! Number of plasma passes - integer :: ipol ! Whether to compute wave polarisation + real(wp_) :: dst ! Integration step size + real(wp_) :: rwmax = 1 ! Normalized maximum radius of beam power + integer :: nray, nrayr, nrayth ! Total/radial/angular number of rays + integer :: igrad = 0 ! Complex eikonal switch + integer :: nstep = 12000 ! Max number of integration steps + integer(kind(step_enum)) :: idst = STEP_TIME ! Choice of the integration variable + integer :: ipass = 1 ! Number of plasma passes + logical :: ipol = .false. ! Whether to compute wave polarisation (from chi, psi) end type ! EC resonant heating & Current Drive parameters type ecrh_cd_parameters - integer :: iwarm ! Choice of power absorption model - integer :: ieccd ! Choice of current drive model - integer :: ilarm ! Larmor-radius expansion order - integer :: imx ! Max iterations for solving the dispersion relation + integer(kind(absorption_enum)) :: iwarm = ABSORP_FULL ! Choice of power absorption model + integer(kind(current_drive_enum)) :: ieccd = CD_NEOCLASSIC ! Choice of current drive model + integer :: ilarm = 5 ! Larmor-radius expansion order + integer :: imx = -20 ! Max iterations for solving the dispersion relation end type ! Output data parameters type output_parameters - integer :: ipec ! Grid spacing for profiles (profili EC) - integer :: nrho ! Number of grid steps for EC profiles + 1 - integer :: istpr ! Subsampling factor for beam cross-section (units 8, 12) - integer :: istpl ! Subsampling factor for outer rays (unit 33) + integer :: ipec = RHO_TOR ! Radial coordinate of the EC profiles + integer :: nrho = 501 ! Number of grid steps for EC profiles + 1 + integer :: istpr = 5 ! Subsampling factor for beam cross-section (units 8, 12) + integer :: istpl = 5 ! Subsampling factor for outer rays (unit 33) end type ! Other parameters @@ -134,7 +218,7 @@ module gray_params end type ! Global variables exposed for gray_core - integer, save :: iequil, iprof, ipol + integer, save :: iequil, iprof integer, save :: iwarm, ilarm, imx, ieccd integer, save :: igrad, idst, ipass integer, save :: istpr0, istpl0 @@ -228,7 +312,7 @@ contains write(strout(17), '("# RFL rwall:",1x,e12.5)') params%misc%rwall ! code parameters - write(strout(18), '("# COD igrad idst ipass ipol:",4(1x,i4))') & + write(strout(18), '("# COD igrad idst ipass ipol:",3(1x,i4),1x,l4)') & params%raytracing%igrad, params%raytracing%idst, & params%raytracing%ipass, params%raytracing%ipol write(strout(19), '("# COD nrayr nrayth nstep rwmax dst:",5(1x,g0.4))') & @@ -259,15 +343,12 @@ contains character(*), intent(in) :: name, value integer(kind(ini_error)) :: err + character(len=:), allocatable :: temp + ! hope for the best... err = ERR_SUCCESS - select case (name) #include "gray_params.inc" - case default - ! unknown parameter - err = ERR_UNKNOWN - end select end function update_parameter @@ -311,7 +392,7 @@ contains integer, intent(out) :: err ! local variables - integer :: u + integer :: u, temp(3) u = get_free_unit() @@ -324,14 +405,20 @@ contains read(u, *) params%raytracing%nrayr, params%raytracing%nrayth, & params%raytracing%rwmax - read(u, *) params%raytracing%igrad, params%raytracing%ipass, & - params%raytracing%ipol + read(u, *) params%raytracing%igrad, params%raytracing%ipass, temp(1) + + ! convert integers to logicals + params%raytracing%ipol = temp(1) > 0 + read(u, *) params%raytracing%dst, params%raytracing%nstep, & params%raytracing%idst read(u, *) params%ecrh_cd%iwarm, params%ecrh_cd%ilarm, params%ecrh_cd%imx read(u, *) params%ecrh_cd%ieccd + ! restrict to 0-3 + params%ecrh_cd%ieccd = min(params%ecrh_cd%ieccd, 3) + read(u, *) params%antenna%alpha, params%antenna%beta read(u, *) params%antenna%power read(u, *) params%antenna%iox, params%antenna%psi, params%antenna%chi @@ -340,8 +427,13 @@ contains read(u, *) params%equilibrium%iequil read(u, *) params%equilibrium%filenm - read(u, *) params%equilibrium%icocos, params%equilibrium%ipsinorm, & - params%equilibrium%idesc, params%equilibrium%ifreefmt + read(u, *) params%equilibrium%icocos, temp(1:3) + + ! convert integers to logicals + params%equilibrium%ipsinorm = temp(1) > 0 + params%equilibrium%idesc = temp(2) > 0 + params%equilibrium%ifreefmt = temp(3) > 0 + read(u, *) params%equilibrium%ixp, params%equilibrium%ssplps, & params%equilibrium%ssplf read(u, *) params%equilibrium%sgnb, params%equilibrium%sgni, & @@ -358,6 +450,9 @@ contains read(u, *) params%output%ipec, params%output%nrho read(u, *) params%output%istpr, params%output%istpl + ! remap to 0,1 (as in irho) + params%output%ipec = mod(params%output%ipec, 2) + close(u) end subroutine read_gray_params @@ -384,7 +479,6 @@ contains mod="gray_params", proc="set_globals") end if - ipol = params%raytracing%ipol igrad = params%raytracing%igrad idst = params%raytracing%idst ipass = params%raytracing%ipass diff --git a/src/gray_params.sh b/src/gray_params.sh deleted file mode 100644 index ae21d72..0000000 --- a/src/gray_params.sh +++ /dev/null @@ -1,30 +0,0 @@ -#!/bin/sh - -# This script generates the very repetitive code for parsing the -# GRAY parameters from a string. The output is embedded -# into gray_params.f90 using a CPP include directive. - -sets='antenna equilibrium profiles raytracing ecrh_cd output misc' - -# shellcheck disable=SC2034 -{ -antenna='alpha beta power psi chi iox ibeam filenm fghz pos w ri phi' -equilibrium='ssplps ssplf factb sgnb sgni ixp iequil icocos ipsinorm idesc ifreefmt filenm' -profiles='psnbnd sspld factne factte iscal irho iprof filenm' -raytracing='rwmax dst nrayr nrayth nstep igrad idst ipass ipol' -ecrh_cd='iwarm ilarm imx ieccd' -output='ipec nrho istpr istpl' -misc='rwall' -} - -deref() { eval "echo \$$1"; } - -for set in $sets; do - for param in $(deref "$set"); do - cat < 0) err = ERR_VALUE -EOF - done -done diff --git a/src/multipass.f90 b/src/multipass.f90 index 46078f9..8a8e056 100644 --- a/src/multipass.f90 +++ b/src/multipass.f90 @@ -165,7 +165,8 @@ contains subroutine initmultipass(i,iox,iroff,yynext,yypnext,yw0,ypw0,stnext,p0ray, & taus,tau1,etau1,cpls,cpl1,lgcpl1,psipv,chipv) ! initialization before pass loop ! arguments - integer, intent(in) :: i, iox ! ipol, mode active on 1st pass + logical, intent(in) :: i ! ipol + integer, intent(in) :: iox ! mode active on 1st pass logical, dimension(:,:), intent(out), pointer :: iroff ! global ray status (F = active, T = inactive) real(wp_), dimension(:), intent(out), pointer :: p0ray,tau1,etau1,cpl1,lgcpl1, & psipv,chipv @@ -173,7 +174,7 @@ contains real(wp_), dimension(:,:,:), intent(out), pointer :: yynext,yypnext ! iroff = .false. ! global ray status (F = active, T = inactive) - if(i.eq.0) call turnoffray(0,1,3-iox,iroff) ! ipol = 0 => stop other mode (iox=1/2 -> stop ib=2/1 at first pass) + if(.not. i) call turnoffray(0,1,3-iox,iroff) ! !ipol => stop other mode (iox=1/2 -> stop ib=2/1 at first pass) yynext = zero ! starting beam coordinates (1) yypnext = zero ! starting beam coordinates (2) yw0 = zero ! temporary beam coordinates (1) diff --git a/src/pec.f90 b/src/pec.f90 index cbaa5ac..91c5e1d 100644 --- a/src/pec.f90 +++ b/src/pec.f90 @@ -26,8 +26,8 @@ contains ! rt_in present: read input grid ! else: build equidistant grid dimension nnd +! ipec=0 rho_tor grid ! ipec=1 rho_pol grid -! ipec=2 rho_tor grid call dealloc_pec allocate(rhop_tab(nnd),rhot_tab(nnd),rtabpsi1(0:nnd),dvol(nnd),darea(nnd), & ratjav(nnd),ratjbv(nnd),ratjplv(nnd)) diff --git a/tests/06-ITER-startup/__init__.py b/tests/06-ITER-startup/__init__.py index 0cb130d..75a812d 100644 --- a/tests/06-ITER-startup/__init__.py +++ b/tests/06-ITER-startup/__init__.py @@ -21,4 +21,4 @@ class Mixed(GrayTest, TestCase): Horizontal polarisation at launch (both X, O mode) ''' reference = get_basedir(__name__) / 'outputs' / 'mixed' - gray_params = {"raytracing.ipol": 1} + gray_params = {"raytracing.ipol": True}