!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 ! D.J. Jeffery, 2020jan06 ! include 'saha_solver.f' ! The Saha solver. !23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 include 'xne_acceleration.f' ! Electron density solution iteration acceleration. include 'partition_function_read.f' ! Reads in the partition function table. !23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 subroutine saha_solver(kkel,kel,imassnumfrac,xel, & & rho,tem,r_dilution, & & iacceleration,idiagnostic,ih_minus_not, & & xne_converged, & & ixne_initial,x_reduced,xx_reduced, & & xne,xne_h,aname,aweight,xne_atom_ratio, & & fraction,iter) use precision_set use contants include 'module_implicit.f' ! Paramenters have to precede their first use. integer, parameter :: nel=38 ! Only to A=38 Sr is allowed at present. integer, parameter :: nion=5 ! Ionization stages, -1,0,1,2,5 are calculated---except for ! ! for those that do not have them, of course. ! ! Stage 5 is treated as having infinite ionization energy ! ! represented by 9999.999 which is just a placeholder. ! ! The stage -1 is calculated for only H: i.e., H-. ! ! The ih_minus_not = 1 turns off the H- ! ! which is useful for comparison to the ! ! the analytic pure H xne_h which does not use stage H-. include 'function_type_declaration.f' ! A file that declares the function types. character, intent(inout) :: aname(nel)*2 ! Output: Atomic name symbols. integer :: i,j,k,kk,L,m,n ! Usually indexes for doloops. However, following ! ! Mihalas-1978-112, j is also ionization stage and k is element sometimes. integer, intent(in) :: iacceleration ! The parameter controling the acceleration of the iteration ! ! for electron density. ! ! 0 for none, 1 for midpoint, 2 for linear, 3 for quadratic. ! ! See discussion of acceleration in ! ! file xne_acceleration.f. ! ! iacceleration=3 recommeded circa 2020jan06. ! https://stackoverflow.com/questions/37723973/fortran-2003-2008-elegant-default-arguments ! According to this site, f95 does NOT allow default values for dummy arguments!!! ! MRC-84 says optional attribute allowed. But how do you tell what to do when ! you do not supply a value? ! The website above has no elegant, universal workaround to suggest. integer, intent(in) :: idiagnostic ! Input: 0 for no diagnostics of the saha solver. ! ! 1 for diagnostics for saha_solver.f ! ! 2 for diagnostics for saha_solver.f ! ! and 'partition_function_read.f' integer, intent(in) :: ih_minus_not ! Input: Set to 0 for H- used. ! ! Set to 1 to turn off H- which is useful for comparison to the ! ! the analytic solution for pure H which does not have ! ! the H- stage. integer :: ih_minus_not_1 ! Equals ih_minus_not, unless ixne_initial=3 in which case it is 1. integer, intent(in) :: imassnumfrac ! Input: 1 for mass fraction, 2 for number fraction. integer, save :: ipartfunc = 0 ! 0 initially causes a read-in of the partition function table. ! ! 1 after the partition function has been read. integer, intent(out) :: iter ! Output: Count of iterations to convergence. integer, intent(in) :: ixne_initial ! Input: For the initial electron density. ! ! 1 the anlytical hydrogen formula in which case the ! ! input xne is ignored. ! ! 2 uses the input xne as the initial value which is a good ! ! idea if you are doing a for series of states where each ! ! succeeding state is close to the last one, and so ! ! the converged last state xne is a good initial xne for ! ! the next state xne. ! ! 3 is for testing with pure H (with no H-) where one supplies ! ! reduced ion density xx_reduced ! ! and an initial reduced electron density ! ! see if on converges to the exact result. ! ! The initial xne is ignored in this case. integer :: jmin,jmax ! Minimum and maximum ionization stages. Which vary with the element. integer, intent(in) :: kkel ! Input: Number of elements to be calculated for. ! ! Except for ixne_initial=3 in which case only H is solved for without H-. integer :: kkel_1 ! Equals kkel, unless ixne_initial=3 in which case it is 1. integer, intent(in) :: kel(kkel) ! Input: Elements by atomic number to be calculated for. real (kind=np), intent(inout) :: aweight(nel) ! Output: Atomic weights from the read-in ! partition function table. ! ! They are rounded off to whole numbers from ! ! terrestrial averages (I think). The exact ! ! isotope weighted atomic weights hard ! ! to know a priori for general astrophysical systems. real (kind=np) :: aweight_1(kkel) ! Atomic weights for the input elements. real (kind=np) :: denom ! A denominator value used briefly. real (kind=np), intent(out) :: fraction(-1:nion,kkel) ! Output: Ionization fractions. real (kind=np), save :: partfunc(-1:nion,nel) ! Initially the partition functions which are currently ! ! just the ground state statistical weights from ! ! the read-in table. But these are then converted ! ! to the ratios of the partition functions since ! ! those are actually used. See Mihalas-1978-113. real (kind=np) :: pair(2,3) ! Used for the linear acceleration of the iteration: ! ! slope=(pair(2,2)-pair(2,1))/(pair(1,2)-pair(1,1)) . real (kind=np) :: phi(-1:nion,kkel) ! Mihalas-1978-113's phi-tilde function. real (kind=np) :: phi_prod(-1:nion,kkel) ! Mihalas-1978-114's phi-tilde function products. real (kind=np), intent(in) :: r_dilution ! Input: If r_dilution < 0, does nothing. ! ! Otherwise tem(used)=tem(input)*W**(1/4), ! ! where w=(1/2)(1-sqrt(1-r_dilution**(-2))) ! ! is the dilution factor (Mihalas-1978-120) ! ! for photospheric approximation for temperature. ! ! r_dilution = radius/radius_photosphere. real (kind=np) :: relerr ! Relative change-in or relative error-in. real (kind=np), intent(in) :: rho ! Input: Matter density in g/cm**3. real (kind=np) :: sumel ! Sum of the element abundances. real (kind=np), intent(in) :: tem ! Input: Temperature in kelvin. real (kind=np) :: tem_use ! Temperature in kelvin used in calculation. ! ! It is tem if r_dilution < 0 and the dilution ! ! temperature if r_dilution > 0. real (kind=np) :: w_dilution ! If r_dilution > 0, it is the diltuion factor. real (kind=np), intent(out) :: xne_atom_ratio ! Output: Ratio (electron density/atom density). real (kind=np), intent(in) :: xne_converged ! Input: The convergence criterion. ! ! The iteration stops when ! ! relative changes in electron density from ! ! xne to xne_out are less/equal to ! ! xne_converged. Probaly xne_converged=1/100 ! ! is good enough for ! ! actual use of the saha solver and not just testing. ! ! We know errors in the converged value to 1 percent ! ! are not the limiting errors. real (kind=np), optional, intent(in) :: x_reduced,xx_reduced ! Input. Values used for ! ! testing in the ixne_initial=3 case. real (kind=np), intent(in) :: xel(kkel) ! Input: Mass/number fraction elements. real (kind=np), save :: xion_en(-1:nion,nel) ! Ionization energies in eV from the read-in ! ! partition function table. real (kind=np) :: xnum(kkel) ! Number density of elements. real (kind=np) :: xnum_total ! Total number density of elements collectively real (kind=np), intent(inout) :: xne ! Input/Output: Electron density, input is inital value ! ! for iteration of ixne_initial=2 and output is the ! ! converged value real (kind=np) :: xne_accel ! The accelerated electron density for the next iteration ! ! of the Saha electron density equation. real (kind=np), intent(out) :: xne_h ! Output: Electron density for the pure hydrogen (without H-) ! ! is used as the initial value for the iteration if xne < 0. ! ! The analytic solution is exact (within numerical error) ! ! for partition functions used. real (kind=np) :: xne_L,xne_q ! The linear and quadratic accelerated xne values for next ! ! Saha equation evaluation. ! ! iacceleration=(0,1,2,3) chooses which of ! ! xne_out,xne_mid,xne_L,xne_q to use. real (kind=np) :: xne_min,xne_max ! The current allowed range for the converged electron density. real (kind=np) :: xne_mid ! The midpoint value for current allowed range for the ! ! converged electron density. real (kind=np) :: xne_out ! The electron density value produced by the Saha electron density equation. real (kind=np) :: xne_j ! The unnormalized electron density contribution of ionization stage j ! ! of element k. ! if(ipartfunc .eq. 0)then ipartfunc=1 call partition_function_read(nel,nion,idiagnostic, & & aname,aweight,xion_en,partfunc) partfunc(-1:nion-1,:)=partfunc(-1:nion-1,:)/partfunc(0:nion,:) ! ! The ratios are actually what is needed (Mihalas-1978-113). if(idiagnostic .eq. 2) then do405 : do k=1,nel ! k is atomic number as for Mihalas-1978-112. jmax=min(k,nion) write(*,fmt=910) aname(k),k,aweight(k), & & xion_en(-1:jmax,k) write(*,fmt=920) partfunc(-1:jmax,k) 910 format(1x,a2,i4,f6.0,18x,7f9.3) 920 format(31x,7f9.3) end do do405 end if end if ! if(iacceleration .lt. 0 .or. iacceleration .gt. 3) then print*,'You have given iacceleration=',iacceleration print*,'Allowed iacceleration values 0,1,2,3' stop 'saha_solver.f' end if ! do406 : do kk=1,kkel ! kk is the calculation index number for the element. aweight_1(kk)=aweight(kel(kk)) end do do406 ! sumel=sum(xel(:)) ! For normalization if that was not set initially. if(imassnumfrac .eq. 1) then ! For mass fractions as input parameters. xnum(:)=rho*(xel(:)/sumel)/(aweight_1(:)*amu) xnum_total=sum(xnum(:)) else ! For number fractions as input parameters. xnum_total=rho/sum((xel(:)/sumel)*aweight_1(:)*amu) xnum(:)=xnum_total*(xel(:)/sumel) end if ! if(r_dilution .le. 0.0_np) then tem_use=tem else w_dilution=0.5_np*(1.0_np-sqrt(1.0_np-r_dilution**(-2))) ! Mihalas-1978-120 tem_use=tem*w_dilution**0.25_np end if do408 : do kk=1,kkel phi(-1:nion-1,kk)=partfunc(-1:nion-1,kel(kk)) & ! Mihalas-1978-113 & *saha_con*tem_use**(-1.5_np) & & *exp(xion_en(-1:nion-1,kel(kk))/(xk_boltzmann_ev*tem_use)) ! print*,'phi assigned, etc.' ! print*,partfunc(0,kel(kk)),saha_con,tem_use ! print*,xion_en(0,kel(kk)), ! & xion_en(0,kel(kk))/(xk_boltzmann_ev*tem_use),phi(0,1) end do do408 ! do410 : do kk=1,kkel if(kel(kk) .gt. 1) then jmin=0 else jmin=-1 end if jmax=min(kel(kk),nion) phi_prod(jmax,kk)=1.0_np do430 : do j=jmax-1,jmin,-1 phi_prod(j,kk)=phi(j,kk)*phi_prod(j+1,kk) ! Mihalas-1978-114 end do do430 end do do410 ! kkel_1=kkel ih_minus_not_1=ih_minus_not if(ixne_initial .eq. 1) then xne_h=func_xne_h_analytic_general(xnum_total*phi(0,1)) & & /phi(0,1) if(idiagnostic .ge. 1) then print*,'xne_h,xnum_total,phi(0,1)' write(*,'(1p,3e15.7)') xne_h,xnum_total,phi(0,1) end if xne=xne_h else if(ixne_initial .eq. 2) then if(xne .le. 0.0_np) stop ! print*,'For ixne_initial=2 you need xne > 0' ! print*,'Your input value is xne= ',xne else xne=x_reduced/phi(0,1) xnum_total=xx_reduced/phi(0,1) kkel_1=1 ! This is a pure hydrogen case no matter what the other parameters are. ih_minus_not_1=1 ! No H- in this case. end if ! if(idiagnostic .ge. 1) then write(*,'(a)') 'The Saha solution for:' write(*,fmt=930) iacceleration,rho,tem_use,xne 930 format('iacceleration=',i5,1p,' rho=',e15.7, & & ' T=',e15.7,' xne_initial=',e15.7) write(*,*) 'kel(:)',kel(:) write(*,*) 'xel(:)',xel(:) write(*,fmt=940) end if iter=0 pair(:,:)=0.0_np do510 : do iter=iter+1 xne_out=0.0_np do520 : do kk=1,kkel_1 if(kel(kk) .gt. 1 .or. ih_minus_not_1 .eq. 1) then jmin=0 else jmin=-1 ! The case where H- is included. end if jmax=min(kel(kk),nion) xne_j=0.0_np do530 : do j=jmin,jmax fraction(j,kk)=xne**(jmax-j)*phi_prod(j,kk) ! Mihalas-1978-114--115 xne_j=xne_j+real(j,kind=np)*fraction(j,kk) ! Mihalas-1978-114--115 ! See MRC-162 for real function. ! print*,fraction(j,kk),xne_j end do do530 denom=sum(fraction(jmin:jmax,kk)) fraction(j,kk)=fraction(j,kk)/denom xne_out=xne_out+xnum(kk)*xne_j/denom end do do520 relerr=(xne_out-xne)/xne if(abs(relerr) .le. xne_converged) then xne_atom_ratio=xne_out/xnum_total if(idiagnostic .ge. 1) then write(*,fmt=942) iter,xne,xne_out,0.0_np,relerr, & & xne_atom_ratio end if xne=xne_out ! Convergence achieved. exit do510 end if pair(:,1)=pair(:,2) pair(:,2)=pair(:,3) pair(1,3)=xne pair(2,3)=xne_out call xne_acceleration(iacceleration,iter, & & pair,xne,xne_out, & & xne_min,xne_max,xne_mid,xne_L,xne_q,xne_accel) if(idiagnostic .ge. 1) then write(*,fmt=942) iter,xne,xne_out,xne_accel,relerr, & & xne_accel/xnum_total if(iter .eq. 100) stop 'The iteration failing iter=100' end if xne=xne_accel ! This must come after the idiagnostic printout. end do do510 ! 940 format(' iter',12x,'xne',8x,'xne_out',6x,'xne_accel',9x,'relerr', & & 5x,'Ionization') 942 format(i5,1p,5e15.7) ! end subroutine saha_solver !23456789a123456789b123456789c123456789d123456789e123456789f123456789g12