diff --git a/pentrc/inputs.f90 b/pentrc/inputs.f90 index d764fffd..053462f9 100644 --- a/pentrc/inputs.f90 +++ b/pentrc/inputs.f90 @@ -201,7 +201,6 @@ subroutine read_kin(file,zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,write_log,format) 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 @@ -209,7 +208,7 @@ subroutine read_kin(file,zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,write_log,format) 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) :: zeff,zpitch,welec,wdian,wdiat real(r8), dimension(0:nkin,2) :: tkev real(r8), dimension(0:nkin,3) :: nu_supplied type(spline_type) :: tmp,tmp_nu @@ -288,22 +287,11 @@ subroutine read_kin(file,zi,zimp,mi,mimp,nfac,tfac,wefac,wpfac,write_log,format) ! manipulation of rotation variables welec(:) = wefac*kin%fs(:,5) ! direct manipulation of omegae - ! note dividing by zero makes a nan that then ruins the whole spline - warning_printed = .false. - do i=0,nkin - if(welec(i) == 0) then - welec(i) = 1e-9 - if(.not. warning_printed)then - print *, " ! WARNING ! Modifying omega_ExB = 0 to 1e-9 to avoid nans" - warning_printed = .true. - end if - end if - end do if(wefac/=1.0 .and. verbose) print('(a55,es10.2e3)'),' -> applying direct manipulation of omegaE by factor ',wefac wdian =-twopi*kin%fs(:,3)*kin%fs1(:,1)/(e*zi*chi1*kin%fs(:,1)) wdiat =-twopi*kin%fs1(:,3)/(e*zi*chi1) - wpefac= (wpfac*(welec+wdian+wdiat) - (wdian+wdiat))/welec - welec = wpefac*welec ! indirect manipulation of rotation + ! wpfac scales welec+wdian+wdiat without dividing by welec, so 0 stays 0 + welec = welec + (wpfac-1)*(welec+wdian+wdiat) kin%fs(:,5)=welec if(wpfac/=1.0 .and. verbose) then diff --git a/regression/test/pentrc_zero_omega_e.f90 b/regression/test/pentrc_zero_omega_e.f90 new file mode 100644 index 00000000..98c91079 --- /dev/null +++ b/regression/test/pentrc_zero_omega_e.f90 @@ -0,0 +1,29 @@ +! Usage: pentrc_zero_omega_e +! read_kin must keep the controls exact: wefac = 0 with wpfac = 1 gives +! omega_E = 0 at every knot, and wpfac = 1 leaves wefac*omega_E unchanged. +program pentrc_zero_omega_e + use params, only: r8 + use inputs, only: read_kin, kin, chi1 + use spline_mod, only: spline_dealloc + implicit none + character(512) :: file + real(r8), allocatable :: welec_full(:) + + call get_command_argument(1, file) + chi1 = 1.0_r8 + + call read_kin(file, 1, 6, 2, 12, 1.0_r8, 1.0_r8, 1.0_r8, 1.0_r8, .false.) + welec_full = kin%fs(:, 5) + if (all(welec_full == 0)) error stop 'fixture has no rotation' + + call spline_dealloc(kin) + call read_kin(file, 1, 6, 2, 12, 1.0_r8, 1.0_r8, 0.0_r8, 1.0_r8, .false.) + print '(a,es12.4)', 'max |omega_E| for wefac = 0: ', maxval(abs(kin%fs(:, 5))) + if (any(kin%fs(:, 5) /= 0)) error stop 'wefac = 0 does not give omega_E = 0' + + call spline_dealloc(kin) + call read_kin(file, 1, 6, 2, 12, 1.0_r8, 1.0_r8, 0.5_r8, 1.0_r8, .false.) + if (any(kin%fs(:, 5) /= 0.5_r8*welec_full)) & + error stop 'wpfac = 1 changes omega_E' + print '(a)', 'PASS: exact omega_E controls' +end program pentrc_zero_omega_e