PENTRC - NEW FEATURE - Optional supplied collision rates (kinetic_format="rates-v1")

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))

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"
 

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,45 @@ 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))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"
+            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 +253,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,20 +264,21 @@ 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) &
             /(sqrt(1.0*mi)*tkev(:,1)**1.5)
         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"

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