From 10e850d381c992a0c47c6d596cc1ad119e9df990 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Mon, 19 Oct 2026 06:14:14 -0400 Subject: [PATCH 1/2] PENTRC - NEW FEATURE - Optional supplied collision rates via kinetic_format="rates-v1" With kinetic_format="rates-v1" in PENT_INPUT, the kinetic_file carries three extra columns (logLambda, nu_i, nu_e in s^-1) that replace the internal collision model. They are interpolated in log on the table so nodes stay positive, and they reuse the pchip nu slopes from the edge-collisionality fix. The default "legacy" path is unchanged. DCON kin_flag passes the same format. Slim replacement for the supplied-rates part of #322 on top of #336. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VMDGyfQKRZQMdk3NJXsKQE --- dcon/dcon.F | 2 +- input/pentrc.in | 1 + pentrc/inputs.f90 | 37 +++++++++++++++++++++++++++++++------ pentrc/pentrc_interface.f90 | 7 ++++--- 4 files changed, 37 insertions(+), 10 deletions(-) diff --git a/dcon/dcon.F b/dcon/dcon.F index df3dcde7e..bd125dc9b 100644 --- a/dcon/dcon.F +++ b/dcon/dcon.F @@ -222,7 +222,7 @@ PROGRAM dcon $ twopi*psio,ro,nn,jac_type,mlow,mhigh,mpert,mthvac) ! manually set the kinetic profiles CALL read_kin(kinetic_file,zi,zimp,mi,mimp,nfac, - $ tfac,wefac,wpfac,indebug) + $ tfac,wefac,wpfac,indebug,kinetic_format) ! manually set the perturbed equilibrium displacements ! use false flat xi and xi' for equal weighting ALLOCATE(psitmp(sq%mx+1),mtmp(mpert),xtmp(sq%mx+1,mpert)) diff --git a/input/pentrc.in b/input/pentrc.in index b1ba6e95c..f6f3f4ba9 100644 --- a/input/pentrc.in +++ b/input/pentrc.in @@ -1,6 +1,7 @@ &PENT_INPUT idconfile = "euler.bin" ! Eigenfunctions output by DCON run kinetic_file = "kinetic.txt" ! Ascii profiles. columns: psi_n, ni(m^-3), ne(m^-3), ti(eV), te(eV), omega_EXB(rad/s) + kinetic_format = "legacy" ! "rates-v1" adds columns logLambda, nu_i(s^-1), nu_e(s^-1) that replace the internal collision model data_dir = "" ! Location of pre-formed fnml matrix used in rlar and clar methods. "default" uses $GPECHOME/pentrc. peq_file = "" ! Ascii displacement. columns: psi_n, m, real(xi^psi), imag(xi^psi), real(xi^alpha), imag(xi^alpha). Default uses "gpec_xclebsch_n#.out" diff --git a/pentrc/inputs.f90 b/pentrc/inputs.f90 index 32d626259..db5f28a02 100644 --- a/pentrc/inputs.f90 +++ b/pentrc/inputs.f90 @@ -156,7 +156,7 @@ end subroutine read_equil !======================================================================= - subroutine read_kin(file,zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,write_log) + subroutine read_kin(file,zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,write_log,format) !----------------------------------------------------------------------- !*DESCRIPTION: ! Read ascii file containing table of kinetic profiles ni, ne, ti, @@ -168,6 +168,10 @@ subroutine read_kin(file,zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,write_log) ! header and/or footer, with the exception that no lines start with ! a number. ! + ! With format="rates-v1" the table has 9 columns, the 6 above plus + ! logLambda, nu_i(s^-1), nu_e(s^-1). These positive supplied values + ! replace the internal collision model (interpolated in log). + ! !*ARGUMENTS: ! file : character(256) (in) ! File path. @@ -185,6 +189,8 @@ subroutine read_kin(file,zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,write_log) ! Scaling of rotation profile, done via manipulation of omegaE ! write_log : bool ! Writes kinetic spline to log file + ! format : character (optional) + ! "legacy" (default) or "rates-v1" ! !----------------------------------------------------------------------- @@ -194,26 +200,43 @@ subroutine read_kin(file,zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,write_log) integer, intent(in) :: zi,zimp,mi,mimp real(r8), intent(in) :: wefac,wpfac,nfac,tfac logical :: write_log + character(*), intent(in), optional :: format ! declare local variables logical :: warning_printed = .false. real(r8), dimension(:,:), allocatable :: table character(32), dimension(:), allocatable :: titles integer, parameter :: nkin = 100 - integer :: i,out_unit, tshape(2) + integer :: i,out_unit, tshape(2), nq + logical :: supplied real(r8) :: psi real(r8), dimension(0:nkin) :: zeff,zpitch,welec,wdian,wdiat,wpefac real(r8), dimension(0:nkin,2) :: tkev + real(r8), dimension(0:nkin,3) :: nu_supplied type(spline_type) :: tmp,tmp_nu call readtable(file,table,titles,verbose,write_log) tshape = shape(table) if(write_log) print *,"Table shaped ",tshape + ! optional supplied logLambda, nu_i, nu_e columns + supplied = .false. + if(present(format)) supplied = (format == "rates-v1") + if(present(format) .and. .not. supplied .and. format /= "legacy") & + error stop "ERROR: read_kin - kinetic_format must be legacy or rates-v1" + nq = 5 + if(supplied)then + if(tshape(2) /= 9) error stop "ERROR: read_kin - rates-v1 needs 9 columns" + if(any(table(:,7:9) <= 0)) & + error stop "ERROR: read_kin - rates-v1 logLambda and nu must be positive" + table(:,7:9) = log(table(:,7:9)) ! interpolate in log to stay positive + nq = 8 + endif + ! temp spline on experimental grid - call spline_alloc(tmp,tshape(1)-1,5) + call spline_alloc(tmp,tshape(1)-1,nq) tmp%xs(0:) = table(1:,1) - do i=1,5 + do i=1,nq tmp%fs(0:,i) = table(1:,i+1) enddo call spline_fit_pchip(tmp) ! monotone between table points @@ -228,6 +251,7 @@ subroutine read_kin(file,zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,write_log) call spline_eval(tmp,psi,0) kin%xs(i) = psi kin%fs(i,1:5) = tmp%f(1:5) + if(supplied) nu_supplied(i,:) = exp(tmp%f(6:8)) enddo if(any(kin%fs(:,1:2) <= 0)) error stop "ERROR: read_kin - densities must be positive" if(any(kin%fs(:,3:4) < 0)) error stop "ERROR: read_kin - temperatures must be non-negative" @@ -238,13 +262,13 @@ subroutine read_kin(file,zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,write_log) ! collisionalities **assumes si units for n,t** ! floor T (keV) at 1 meV and logLambda at 1 so cold (T -> 0) edges give large positive nu tkev = max(kin%fs(:,3:4)/1.602e-16, 1.0e-6_r8) - if(any(kin%fs(:,3:4)/1.602e-16 < 1.0e-6_r8)) & + if(any(kin%fs(:,3:4)/1.602e-16 < 1.0e-6_r8) .and. .not. supplied) & print *, " ! WARNING ! Flooring T at 1e-3 eV in collision frequencies" zeff = zimp-(kin%fs(:,1)/kin%fs(:,2))*zi*(zimp-zi) zpitch = 1.0+(1.0+mimp)/(2.0*mimp)*zimp*(zeff-1.0)/(zimp-zeff) kin%fs(:,6) = 17.3-0.5*log(kin%fs(:,2)/1.0e20) & +1.5*log(tkev(:,2)) - if(any(kin%fs(:,6) < 1)) & + if(any(kin%fs(:,6) < 1) .and. .not. supplied) & print *, " ! WARNING ! Flooring logLambda at 1 in collision frequencies" kin%fs(:,6) = max(kin%fs(:,6), 1.0_r8) kin%fs(:,7) = (zpitch/3.5e17)*kin%fs(:,1)*kin%fs(:,6) & @@ -252,6 +276,7 @@ subroutine read_kin(file,zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,write_log) kin%fs(:,8) = (zpitch/3.5e17)*kin%fs(:,2)*kin%fs(:,6) & /(sqrt(me/mp)*tkev(:,2)**1.5) kin%fs(:,9) = zeff + if(supplied) kin%fs(:,6:8) = nu_supplied call spline_fit(kin,"extrap") if(write_log) print *,"Formed kin spline" diff --git a/pentrc/pentrc_interface.f90 b/pentrc/pentrc_interface.f90 index 6814a986c..2bdbe6de0 100644 --- a/pentrc/pentrc_interface.f90 +++ b/pentrc/pentrc_interface.f90 @@ -143,13 +143,14 @@ module pentrc_interface pmodb_file ="none", & data_dir ="." character(32) :: & + kinetic_format = "legacy",& nutype = "harmonic",& f0type = "maxwellian",& jac_in = "default",& moment = "pressure" ! namelists - namelist/pent_input/kinetic_file, gpec_file, peq_file, pmodb_file, idconfile, & + namelist/pent_input/kinetic_file, kinetic_format, gpec_file, peq_file, pmodb_file, idconfile, & data_dir, zi, zimp, mi, mimp, nl, electron, nutype, f0type, & jac_in, jsurf_in, tmag_in, power_bin, power_bpin, power_rin, power_rcin @@ -238,7 +239,7 @@ subroutine initialize_pentrc(op_kin,op_deq,op_peq) ! read (perturbed) equilibrium inputs if(indebug .and. in_kin) print *," read_kin args: ",trim(kinetic_file),zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,indebug if(in_deq) call read_equil(idconfile) - if(in_kin) call read_kin(kinetic_file,zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,indebug) + if(in_kin) call read_kin(kinetic_file,zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,indebug,kinetic_format) if(in_peq)then if(peq_file=="")then write(nstr,'(i3)')nn @@ -297,7 +298,7 @@ subroutine get_pentrc(get_nl,get_zi,get_mi,get_wdfac,get_divxfac,& print *,"read_kin args:" print *," ",trim(kinetic_file),zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,tdebug endif - call read_kin(kinetic_file,zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,indebug) + call read_kin(kinetic_file,zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,indebug,kinetic_format) ! transfer torque function variables to external program get_nl = nl From 8b6f0d61beb7c5bddbf708e64c9bc9633b13dd1a Mon Sep 17 00:00:00 2001 From: logan-nc Date: Mon, 19 Oct 2026 23:39:01 -0400 Subject: [PATCH 2/2] PENTRC - MINOR - Only test kinetic_format when it is present Fortran does not guarantee short-circuit evaluation, so the combined present(format) .and. format /= "legacy" test could read an absent argument. Every caller passes it today; this nests the check. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VMDGyfQKRZQMdk3NJXsKQE --- pentrc/inputs.f90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/pentrc/inputs.f90 b/pentrc/inputs.f90 index db5f28a02..4d4f2b7ea 100644 --- a/pentrc/inputs.f90 +++ b/pentrc/inputs.f90 @@ -221,9 +221,11 @@ subroutine read_kin(file,zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,write_log,format) ! optional supplied logLambda, nu_i, nu_e columns supplied = .false. - if(present(format)) supplied = (format == "rates-v1") - if(present(format) .and. .not. supplied .and. format /= "legacy") & - error stop "ERROR: read_kin - kinetic_format must be legacy or rates-v1" + if(present(format))then + supplied = (format == "rates-v1") + if(.not. supplied .and. format /= "legacy") & + error stop "ERROR: read_kin - kinetic_format must be legacy or rates-v1" + endif nq = 5 if(supplied)then if(tshape(2) /= 9) error stop "ERROR: read_kin - rates-v1 needs 9 columns"