!23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 ! D.J. Jeffery, 2020jan06 ! !23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 ! Includes standard modules: precision,numerical,others eventually include 'module_standard.f' include 'saha_solver.f' ! The Saha solver. ! include 'xne_acceleration.f' ! Electron density solution iteration acceleration. ! This subroutine is included in saha_solver.f. include 'xne_h_analytic.f' ! The analytic H subprograms. ! ! H- is not included in this analytic solution. !23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 ! program test use precision_set use contants ! use ieee_exceptions ! My compiler does not support it. MRC-221. ! use ieee_arithmetic ! My compiler does not support it. MRC-221. include 'test_declarations.f' ! ! Execution lines begin ! ! 1) Tests of precision and range of types and some variable values. print* print*,'1) Tests of precision and range of types' print*,' and some variable values.' print*,'iprecision_np' print*,iprecision_np print*,'iminexponent_np,imaxexponent_np' print*,iminexponent_np,imaxexponent_np print*,'iprecision_ns' print*,iprecision_ns print*,'iminexponent_ns,imaxexponent_ns' print*,iminexponent_ns,imaxexponent_ns print*,'phi_0_golden_ratio,phi_1_golden_ratio' print*,phi_0_golden_ratio,phi_1_golden_ratio print*,'e_rydberg_ev,e_rydberg_ev_reduced' print*,e_rydberg_ev,e_rydberg_ev_reduced ! print*,'saha_con,saha_con_inverse' print*,saha_con,saha_con_inverse print*,' saha_con = constant used in ' print*,' Saha ionization stage ratios' print*,' (Mihalas-1978-113: his value is 2.07e-16 cm**3.)' tmp=2.0_np*(pi_sqroot*alpha_fine_structure*xm_e & & *clight/h_planck)**3 print*,'tmp,saha_con_inverse_natural,' print*,'tmp,saha_con_inverse_natural,temperature_rydberg' print*,tmp,saha_con_inverse_natural,temperature_rydberg ! ! flag=ieee_is_nan(1.0_np) ! Compiler does not support it. MRC-221. ! if(flag .eqv. .true.) stop ! if(flag .eqv. .false.) stop 'flag stop' ! This works. if(flag .eqv. .false.) print*,'Logical variables work.' ! ! 2) Tests for the analytic H ionization or Saha function. print* print*,'2) Tests for the analytic H ionization' print*,' or Saha function.' do410 : do i=-10,10 xni_reduced=10.0_np**(i) xne_e=func_xne_h_analytic_exact(xni_reduced) ! Exact formula numerically good. xne_e2=func_xne_h_analytic_exact2(xni_reduced)! Exact formula numerical bad. xne_small=func_xne_h_analytic_series_small(xni_reduced) xne_large=func_xne_h_analytic_series_large(xni_reduced) ! Series result to 3rd order ! ! for an altnerating series ! ! where the error is less than ! ! last term neglected (Arf-294). relerr=(xne_e-xne_small)/xne_small tmp=(xne_e-xne_large)/xne_large tmp2=(xne_e2-xne_small)/xne_small write(*,'(i5,8e15.7)') i,xni_reduced,xne_e,relerr,tmp, & & xne_e2,tmp2 end do do410 ! stop ! ! 3) Tests for convergence: iteration,mipoint,linear,quadratic' print* print*,'3) Tests for convergence:' print*,' iteration,mipoint,linear,quadratic' steps=3.0_np steps_factor=10.0_np**(1.0_np/steps) xne_initial_min=10.0_np**(-3) xne_initial_max=10.0_np**3 isteps=log10(xne_initial_max/xne_initial_min) & & /log10(steps_factor) + 1.0_np ! The math for isteps comes from B=Af**isteps xni_reduced=1.0_np ! Exact hydrogenic Saha equation solution for this is (golden ratio)-1 ! ! and this our converged solution. do412 : do iacceleration_index=0,3 ! Go through 0--3 acceleration cases. print* write(*,'(a,i2)'),'Acceleration case ',iacceleration_index xne_initial=xne_initial_min xne_iter_prct=0.0_np ! Average number of iterations to xne_prct for each acceleration case. xne_iter_mpre=0.0_np ! Average number of iterations to xne_mpre for each acceleration case. do420 : do i=1,isteps write(*,'(a,i2,a,i2)') & & ' Iteration ',i, & & ' for iacceleration=',iacceleration_index if(iacceleration_index .ge. 3) write(*,fmt=910) 910 format(' Iter',12x,'xne',8x,'xne_out',8x,'xne_mid', & & 10x,'xne_L',10x,'xne_q',6x,'xne_accel',9x,'relerr') xne=xne_initial pair(:,:)=0.0_np ! Initiate all pairs to zero. iter_prct=0 iter_mpre=0 iter=0 do430 : do iter=iter+1 xne_out=func_xne_h_analytic_equation(xne,xni_reduced) ! The H Saha equation, not its solution. ! ! It is used here as iteration equation. relerr=(xne_out-xne)/xne if(abs(relerr) .le. xne_mpre) then xne=xne_out ! Convergence achieved exit do430 end if pair(:,1)=pair(:,2) pair(:,2)=pair(:,3) pair(1,3)=xne pair(2,3)=xne_out ! print*,'Before calling xne_acceleration' call xne_acceleration(iacceleration_index, & iter,pair,xne,xne_out, & & xne_min,xne_max,xne_mid,xne_L,xne_q,xne_accel) if(iter_prct .eq. 0 .and. & & abs(xne_accel-phi_1_golden_ratio)/phi_1_golden_ratio & & .le. xne_prct) iter_prct=iter if(iter_mpre .eq. 0 .and. & & abs(xne_accel-phi_1_golden_ratio)/phi_1_golden_ratio & & .le. xne_mpre) iter_mpre=iter if(iacceleration_index .ge. 3) & & write(*,'(1p,i5,8e15.7)') iter,xne,xne_out,xne_mid, & & xne_L,xne_q,xne_accel,relerr xne=xne_accel ! For the next evaluation of the electron density equation. if(iter .eq. 50) then print*,'iacceleration case ',iacceleration_index stop 'Exceeded allowed iterations for test cases.' end if end do do430 write(*,920) iacceleration_index,iter,iter_prct,iter_mpre 920 format(2x,'iacceleration=',i2,' iteration count=',i3, & & ' iter_prct=',i3,' iter_mpre=',i3) xne_iter_prct=xne_iter_prct+real(iter_prct,kind=np) xne_iter_mpre=xne_iter_mpre+real(iter_mpre,kind=np) xne_initial=xne_initial*steps_factor end do do420 write(*,930) iacceleration_index, & & xne_iter_prct/real(isteps,kind=np), & & xne_iter_mpre/real(isteps,kind=np) 930 format(' Results case iacceleration=',i2, & & ': Avg iter to prct=',f6.3, & & ': Avg iter to mpre=',f6.3) end do do412 stop ! ! 4) Use saha_solver just to read-in and print out the partition functions, ! both from partition_function_read.f and saha_solver_f to check agreement. print*, print*,'4) Use saha_solver just to read-in' print*,' and print out the partition functions.' idiagnostic=3 include 'saha_solver_call.f' ! The Saha solver. ! ! 5) Use saha_solver for pure H (without H-). print* print*,'5) Use saha_solver for pure H (without H-).' tem=7.e+3_np kkel=1 kel(1)=1 xel(1)=1.0_np idiagnostic=2 ih_minus_not=1 ! Means do not use H-. xne_converged=xne_mpre ixne_initial=1 idefault=0 include 'saha_solver_call.f' ! The Saha solver. ! xne_h,xnum_total,phi(0,1) ! 4.9524433E+10 6.0221408E+10 4.3613598E-12 ! The Saha solution for: ! iacceleration= 3 rho= 1.0000000E-13 T= 7.0000000E+03 xne_initial= 4.9524433E+10 ! kel(:) 1 ! xel(:) 1.0000000000000000000 ! iter xne xne_out xne_accel relerr Ionization ! 1 4.9524433E+10 4.9524433E+10 0.0000000E+00 -3.0088504E-19 8.2237256E-01 ! xne_solution=xne ! Use the output xne of the last test and the huge ! ! range for analytic hydrogen test above. write(*,"('The solution xne=',1p,e15.7)") xne_solution xne_tmp=xne_solution*xne_initial_min ! Use the output xne of the last test and the huge ! ! range for analytic hydrogen test above. do510 : do i=1,isteps idiagnostic=1 ih_minus_not=1 ixne_initial=2 print* xne=xne_tmp write(*,'(a,i2,a,i2,1p,2(a,e15.7))') 'Calculation=',i, & & ' with iacceleration=',iacceleration, & & ' and solution xne=',xne_solution, & & ' for initial xne=',xne idefault=0 include 'saha_solver_call.f' ! The Saha solver. xne_tmp=xne_tmp*steps_factor end do do510 ! ! 6) Solve for Carbon ionization fractions. ! http://cococubed.asu.edu/code_pages/eos_ionize.shtml for calculated results. ! print* print*,'6) Solve for Carbon ionization fractions.' kkel=1 kel(1)=6 ! Pure carbon. xel(1)=1.0_np steps=3.0_np steps_factor=10.0_np**(1.0_np/steps) rho=1.0e-4_np tem_min=1.0e+3_np tem_max=1.0e+7_np isteps=log10(tem_max/tem_min) & & /log10(steps_factor) + 1.0_np tem=tem_min ! do610 : do i=1,isteps idefault=0 idiagnostic=2 ixne_initial=2 include 'saha_solver_call.f' ! The Saha solver. write(*,942) i,tem,fraction(0:6,1) tem=tem*steps_factor end do do610 942 format(i5,1p,8e15.7) ! end program test !23456789a123456789b123456789c123456789d123456789e123456789f123456789g12