!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 ! D.J. Jeffery, 2020jan06 ! An acceleration subroutine for solving the Saha electron density equation. ! include 'xne_acceleration.f' ! xne = number density of electrons in my jargon. ! x = reduced number density of electrons in my jargon. ! Since this code fragment works for either, I use "x" inside, but use "xne" ! as part of the code fragment name. I've compromised. !23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 include 'xne_quadratic.f' !23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 ! subroutine xne_acceleration(iacceleration,iter, & & pair,x,x_out, & & x_min,x_max,x_mid,x_L,x_q,x_accel) use precision_set use contants include 'module_implicit.f' ! implicit none integer, intent(in) :: iacceleration ! 0 for none, 1 for midpoint, 2 for linear, 3 for quadratic. integer, intent(in) :: iter ! The current iteration count. real (kind=np), intent(in) :: pair(2,3) ! The last three (input,output) pairs for the ! ! the Saha electron density equation evaluation. real (kind=np) :: slope ! The slope for the linear acceleration. real (kind=np), intent(in) :: x ! The input for last the Saha electron density equation evaluation. real (kind=np), intent(in) :: x_out ! The output for the last Saha electron density equation evaluation. real (kind=np), intent(out) :: x_accel ! The accelerated electron density input for the next ! ! Saha electron density equation evaluation. real (kind=np), intent(inout) :: x_min,x_max ! Range limits on the true x ! ! whatever that is. ! ! These can always be found ! ! since the Saha electron density equation ! ! is strictly decreasing (I think), ! ! except perhaps for x input = infinity. real (kind=np), intent(out) :: x_mid ! midpoint or average acceleration value ! ! which may not be any acceleration at all ! ! but is a safe fall back choice. real (kind=np), intent(out) :: x_L ! Linear fit acceleration value. real (kind=np), intent(out) :: x_q ! Quadratic fit acceleration value ! ! print*,'In xne_acceleration' x_min=min(x,x_out) x_max=max(x,x_out) x_mid=0.5_np*(x_min+x_max) ! The allowed range mipoint. ! print*,'Before no-acceleration and midpoint acceleration.' ! print*,'iacceleration =',iacceleration ! print*,'x_min,x_max,x_mid,x_accel' ! print*,x_min,x_max,x_mid,x_accel if(iacceleration .eq. 0) then ! The straight iteration case which x_accel=x_out ! I believe will always converge ! ! but very slowly. ! ! But it's actually only guaranteed ! ! not to diverege, I think. ! From the analytic H case: ! Results case iacceleration= 0: Avg iter to prct= 6.053: Avg iter to mpre=41.632 ! So 6.053 iterations to converge to 1 %, and 41.632 iterations to converge to machine precision. ! else ! print*,'Before the x_accel assignment.' x_accel=x_mid ! Which is probably no acceleration, but it ! ! is guaranteed to converge, and so ! ! but its a fail-safe if the linear or ! ! quadratic acceleration gives a NaN or is out ! ! of the known range. ! From the analytic H case: ! Results case iacceleration= 1: Avg iter to prct= 6.842: Avg iter to mpre=36.158 ! So 6.842 iterations to converge to 1 %, and 36.158 iterations to converge to machine precision. ! Except that it's guaranteed to coverge, it is probably no better than straight iteration. ! end if ! ! print*,'Before linear and quadratic acceleration.' if(iter .ge. 2) then slope=(pair(2,2)-pair(2,1))/(pair(1,2)-pair(1,1)) x_L=(pair(2,1)-slope*pair(1,1))/(1.0_np-slope) if(iacceleration .ge. 2 .and. & & (x_L .ge. x_min .and. x_L .le. x_max)) x_accel=x_L ! From the analytic H case: ! Results case iacceleration= 2: Avg iter to prct= 6.579: Avg iter to mpre=12.737 ! So 6.579 iterations to converge to 1 %, and 12.737 iterations to converge to machine precision. ! For convergence to 1 %, it's no better than straight iteration and midpoint acceleration. ! But it does to converge to machine precision much faster. ! The x_mid is the fail-safe value if x_L has become a NaN which can happen ! as machine precision is approached. ! My compiler does not allow direct tests for NaN's, but they cannot be true for ! if statements. ! else x_L=0.0_np slope=0.0_np end if if(iter .ge. 3) then ! print*,'Before calling xne_quadratic' call xne_quadratic(pair,x_q) if(iacceleration .ge. 3 .and. & & (x_q .ge. x_min .and. x_q .le. x_max)) x_accel=x_q ! From the analytic H case: ! Results case iacceleration= 3: Avg iter to prct= 4.105: Avg iter to mpre= 8.211 ! So 4.105 iterations to converge to 1 %, and 8.211 iterations to converge to machine precision. ! This is a significant improvement over the other cases. ! I think that for large numbers of elements and ionization stages, the extra time overhead ! of quadratic acceleration is worthwhile. Also, it's elegant. ! If the iteration starts with a good initial electron density, covergence to of order 1 % ! might be done by the iteration 2 and the quadratic acceleration will never be called. ! Or it might be called just once allowing convergence on the iteration 4. ! I think it's worth the cpu time overhead. My guess is that quadratic iteration ! costs about the cpu time of including 2 or 3 elements. So if one is including 20 ! or more, it's probably worth it. ! The x_mid or x_L values are fail-safe values if x_q has become a NaN which can happen ! as machine precision is approached. ! My compiler does not allow direct tests for NaN's, but they cannot be true for ! if statements. ! else x_q=0.0_np end if ! end subroutine xne_acceleration !23456789a123456789b123456789c123456789d123456789e123456789f123456789g12