diff --git a/equil/direct.f b/equil/direct.f index 69703590..9c423387 100644 --- a/equil/direct.f +++ b/equil/direct.f @@ -13,6 +13,25 @@ c 5. direct_fl_der. c 6. direct_refine. c 7. direct_output. +c 8. find_fl_surface. +c 9. find_fl_surface2. +c 10. direct_xpt_y_out. +c 11. direct_local_xpoint. +c 12. direct_saddle_angle. +c 13. direct_psisaddle. +c 14. direct_xpoint. +c 15. direct_saddle_coords. +c 16. direct_saddle_coords_inv. +c 17. int_validity_check. +c 18. direct_analytic_ints. +c 19. interp_through_zero. +c 20. patch_two_a. +c 21. patch_two_b. +c 22. patch_two_c. +c 23. yinterim. +c 24. analytic_y_out. +c 25. direct_Blocal. +c 26. plot_xpt_convergence. c----------------------------------------------------------------------- c subprogram 0. direct_mod. c module declarations. @@ -27,16 +46,32 @@ MODULE direct_mod INTEGER, PRIVATE :: istep REAL(r8) :: rmin,rmax,zmin,zmax,rs1,rs2 + REAL(r8), DIMENSION(1:2) :: xpt_etas, rxs, zxs, xpt_c11s + REAL(r8), DIMENSION(1:2) :: xpt_gammas, xpt_varthetas + REAL(r8), DIMENSION(1:2) :: xpt_gammas2, xpt_varthetas2 + LOGICAL, DIMENSION(1:2) :: outside_sep + REAL(r8), DIMENSION(2,2) :: xpt_brackets TYPE(bicube_type) :: psi_in LOGICAL :: direct_infinite_loop_flag INTEGER :: direct_infinite_loop_count = 2000 + INTEGER :: num_xpts + + LOGICAL :: debug_xpt,plot_convergence + INTEGER :: btsc1=10 ! used in plot_convergence + INTEGER :: btsc2=30 ! used in plot_convergence + INTEGER :: psfc1=10 ! used in plot_convergence + INTEGER :: psfc2=30 ! used in plot_convergence + LOGICAL :: onecase=.FALSE. ! set true to print local x-point ints TYPE :: direct_bfield_type REAL(r8) :: psi,psir,psiz,psirz,psirr,psizz,f,f1,p,p1 REAL(r8) :: br,bz,brr,brz,bzr,bzz END TYPE direct_bfield_type - REAL(r8) :: etol=1e-8 + REAL(r8) :: etol=1e-8, dqdeps_tol=40, BpBt_tol=0.01, r_tol=1e-5 + REAL(r8) :: xpt_tol=0.2 + INTEGER :: nstepd=20000 + INTEGER :: nstep2=2048 CONTAINS c----------------------------------------------------------------------- @@ -49,23 +84,40 @@ MODULE direct_mod SUBROUTINE direct_run INTEGER :: ir,iz,itheta,ipsi - INTEGER, PARAMETER :: nstep=2048 + INTEGER :: len_y_out1,len_y_out + INTEGER, DIMENSION(2) :: x_starts REAL(r8) :: f0fac,f0,ffac,rfac,eta,r,jacfac,w11,w12,delpsi,q - REAL(r8), DIMENSION(0:nstep,0:4) :: y_out + REAL(r8), DIMENSION(0:(4*nstepd+2*nstep2+6),0:4) :: y_out + REAL(r8), DIMENSION(0:(nstepd+2),0:4) :: y_out1 REAL(r8), DIMENSION(2, mpsi+1) :: xdx REAL(r8), DIMENSION(3,3) :: v + REAL(r8), DIMENSION(1:4) :: difs + REAL(r8), DIMENSION(1:17) :: outmat2f + + LOGICAL :: debug REAL(r8) :: xm,dx,rholow,rhohigh + REAL(r8) :: eval_BpOnBt,eval_dy1,eval_xpt_tol TYPE(direct_bfield_type) :: bf TYPE(spline_type) :: ff + + plot_convergence=.False. + debug=.FALSE. !dumps all spline info to csv + debug_xpt=.TRUE. !verbose x-point outputs + + xpt_etas=0.0 + xpt_brackets=0.0 + num_xpts=0 c----------------------------------------------------------------------- c warning. c----------------------------------------------------------------------- direct_flag=.TRUE. - IF(psihigh >= 1-1e-6)WRITE(*,'(1x,a,es10.3,a)') + IF(psihigh >= 1-1e-6)WRITE(*,'(1x,a,e17.10,a)') $ "Warning: direct equilibrium with psihigh =",psihigh, $ " could hang on separatrix." direct_infinite_loop_flag = .FALSE. + IF(nstepd>30000)CALL program_stop( + $ "Random pointer errors may appear for large arrays..") c----------------------------------------------------------------------- c fit input to cubic splines and diagnose. c----------------------------------------------------------------------- @@ -125,13 +177,59 @@ SUBROUTINE direct_run IF(bin_fl)CALL bin_open(bin_2d_unit,"flint.bin","UNKNOWN", $ "REWIND","none") c----------------------------------------------------------------------- -c start loop over flux surfaces and integrate over field line. +c check for and initialise x-points if present. +c----------------------------------------------------------------------- + !Integrated around psi=1 anticlockwise, auto exits calling + !direct_fl_int if BpBt_tol or xpt_tol exceeded. + CALL direct_fl_int(1.d0-1.d-14,zero,twopi,y_out1,len_y_out1 + $ ,eval_BpOnBt,eval_dy1,eval_xpt_tol,.TRUE.,bf, .FALSE.) + IF(num_xpts>0 .AND. xpt_etas(1) direct_infinite_loop_count) THEN - direct_infinite_loop_flag = .TRUE. - CALL program_stop("Took too many steps to find flux surf.") - ENDIF - ENDDO + CALL find_fl_surface(psifac,eta1,r,z) + CALL direct_get_bfield(r,z,bf,1) psi0=bf%psi c----------------------------------------------------------------------- c initialize variables. c----------------------------------------------------------------------- istep=0 - eta=0 + eta=eta1 deta=twopi/mtheta y=0 y(2)=SQRT((r-ro)**2+(z-zo)**2) @@ -491,13 +586,13 @@ SUBROUTINE direct_fl_int(ipsi,y_out,bf) atol=etol*y(2) iwork=0 rwork=0 - rwork(1)=twopi + rwork(1)=eta2 rwork(11)=0 c----------------------------------------------------------------------- c write header. c----------------------------------------------------------------------- IF(out_fl)THEN - WRITE(out_2d_unit,10)ipsi,psifac + WRITE(out_2d_unit,51)psifac WRITE(out_2d_unit,20) ENDIF c----------------------------------------------------------------------- @@ -512,32 +607,102 @@ SUBROUTINE direct_fl_int(ipsi,y_out,bf) y_out(istep,:)=(/eta,y/) err=(bf%psi-psi0)/bf%psi c----------------------------------------------------------------------- +c compute ratio of Bp to Bt, as well as divergent integrand related +c to q (dy1) +c----------------------------------------------------------------------- + bp=SQRT(bf%br**2+bf%bz**2) + bt=bf%f/r + dy1=rfac/(bf%bz*COS(eta)-bf%br*SIN(eta)) +c----------------------------------------------------------------------- c compute and print output for each step. c----------------------------------------------------------------------- IF(out_fl)WRITE(out_2d_unit,30) $ istep,eta,rwork(11),y(1:2),r,z,bf%psi,err IF(bin_fl)WRITE(bin_2d_unit) - $ REAL(eta,4),REAL(y(1:2),4),REAL(r,4),REAL(z,4),REAL(err,4) + $ REAL(eta,4),REAL(rwork(11),4),REAL(y(1:4),4),REAL(r,4), + $ REAL(z,4),REAL(psifac,4),REAL(err,4) +c----------------------------------------------------------------------- +c xpoint stopping conditions. +c----------------------------------------------------------------------- + IF((bp/bt)dqdeps_tol.AND.psifac>0.9 + $ .AND. (.NOT. simp))THEN + IF(probe_xpt)THEN + !looking for x-points, exit because you've found one + !will automatically call int_validity_check + EXIT + ELSEIF(cooloff)THEN + !cooloff active, do nothing unless you've gone far enough + !to reset it. + IF(ABS(eta-cooleta)>0.2*pi)THEN + cooloff=.FALSE. + ENDIF + ELSE + CALL int_validity_check(r,z,rfac,eta,valid,eval_xpt_tol) + IF(valid)THEN + !its a valid time to use analytic formulas, exit + !numerical integration + EXIT + ELSE!not a valid time to use formulas, power through + !and hope you don't hit one of the accuracy limits + cooloff = .TRUE. + cooleta = eta + ENDIF + ENDIF + ENDIF c----------------------------------------------------------------------- -c advance differential equations. +c regular stopping conditions. c----------------------------------------------------------------------- - IF(eta >= twopi .OR. istep >= nstep .OR. istate < 0 + IF((eta2 > eta1 .AND. eta >= eta2) .OR. + $ (eta2 < eta1 .AND. eta <= eta2) + $ .OR. istep >= nstepd .OR. istate < 0 $ .OR. ABS(err) >= 1)EXIT +c----------------------------------------------------------------------- +c advance differential equations. +c----------------------------------------------------------------------- istep=istep+1 - CALL lsode(direct_fl_der,neq,y,eta,twopi,itol,rtol,atol, + CALL lsode(direct_fl_der,neq,y,eta,eta2,itol,rtol,atol, $ itask,istate,iopt,rwork,lrw,iwork,liw,jac,mf) ENDDO + len_y_out = istep IF(out_fl)WRITE(out_2d_unit,20) IF(bin_fl)WRITE(bin_2d_unit) c----------------------------------------------------------------------- -c abort if istep > nstep. +c abort if istep > nstepd. +c----------------------------------------------------------------------- + IF((bp/bt)dqdeps_tol.AND.psifac>0.9 + $ .AND. (.NOT. simp))THEN + IF(probe_xpt)THEN + rx=r + zx=z + CALL direct_xpoint(rx,zx,num_xpts+1,new_xpt) + IF(new_xpt)THEN + num_xpts=num_xpts+1 + ENDIF + ENDIF c----------------------------------------------------------------------- - IF(eta < twopi)THEN - WRITE(message,40)"direct_int: istep = nstep = ",nstep, - $ " at eta = ",eta,", ipsi = ",ipsi - CALL program_stop(message) +c get tolerance terms. +c----------------------------------------------------------------------- + CALL int_validity_check(r,z,rfac,eta,valid,eval_xpt_tol) + eval_dy1=ABS(dy1) + eval_BpOnBt=(bp/bt) + RETURN + ELSEIF(eta < eta2 .AND. cooloff)THEN + WRITE(message,61)istep,nstepd + WRITE(message2,11)eta,eta1,eta2 + PRINT "(A)", message + PRINT "(3A)","Field line integrator failed before the", + $"analytic x-point formulas became valid. Increase nstepd to ", + $"keep accuracy, or decrease etol/reduce psihigh to retain speed." + CALL program_stop(message2) ENDIF c----------------------------------------------------------------------- +c if .NOT.(eval_dy1==eval_BpOnBt==-1.0), and early exit, we know the +c separatrix integral died because of min_BpOnBt,max_dy1 criteria. +c----------------------------------------------------------------------- + eval_dy1=-1.d0 + eval_BpOnBt=-1.d0 + eval_xpt_tol=-1.d0 +c----------------------------------------------------------------------- c terminate. c----------------------------------------------------------------------- RETURN @@ -706,4 +871,2103 @@ SUBROUTINE direct_output c----------------------------------------------------------------------- RETURN END SUBROUTINE direct_output +c----------------------------------------------------------------------- +c subprogram 8. find_fl_surface. +c finds r,z given psi, eta. Uses a Newton method, similar to +c direct_position. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE find_fl_surface(psifac,eta,r,z) + + REAL(r8), INTENT(IN) :: psifac,eta + REAL(r8), INTENT(OUT) :: r,z + + TYPE(direct_bfield_type) :: bf + REAL(r8), PARAMETER :: eps=1e-13 + INTEGER :: ir + REAL(r8) :: cosfac,sinfac,radius,dradius,dfdradius,psi0 +c----------------------------------------------------------------------- +c find flux surface. +c----------------------------------------------------------------------- + cosfac=COS(eta) + sinfac=SIN(eta) + + psi0=psio*(1-psifac) + radius=SQRT(psifac)*(rs2-ro) + r=ro+cosfac*radius + z=zo+sinfac*radius + ir = 0 + DO + CALL direct_get_bfield(r,z,bf,1) + dfdradius = -bf%psir*cosfac-bf%psiz*sinfac + dradius = -(psi0-bf%psi)/dfdradius + + radius=radius+dradius + r=ro+cosfac*radius + z=zo+sinfac*radius + IF(ABS(dradius) <= eps*r)EXIT + + ir = ir+1 + IF (ir > direct_infinite_loop_count) THEN + PRINT"(A)","find_fl_surface failed, using find_fl_surface2" + CALL find_fl_surface2(psifac,eta,radius,r,z) + EXIT + ENDIF + ENDDO +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE find_fl_surface +c----------------------------------------------------------------------- +c subprogram 9. find_fl_surface2. +c finds r,z given psi, eta. Uses a bisection method instead of the +c Newton method used in find_fl_surface, motivated by special cases +c where the Newton method fails to converge eg. near x-points. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE find_fl_surface2(psifac,eta,radius1,r,z) + + REAL(r8), INTENT(IN) :: psifac,eta,radius1 + REAL(r8), INTENT(OUT) :: r,z + + TYPE(direct_bfield_type) :: bf + REAL(r8), PARAMETER :: eps=1e-13 + INTEGER :: ir + REAL(r8) :: cosfac,sinfac,dradius1,dfdradius,psi0,pos + REAL(r8) :: rbrack_l,rbrack_h,dradius2,dradius3 + + REAL(r8), DIMENSION(1:3) :: radii + + pos=one + rbrack_l=1e-1 + rbrack_h=1e-5 +c----------------------------------------------------------------------- +c find flux surface. +c----------------------------------------------------------------------- + cosfac=COS(eta) + sinfac=SIN(eta) + + psi0=psio*(1-psifac) +c----------------------------------------------------------------------- +c finding lower bound within chosen surface. +c dradius1 should always be greater than 0. +c----------------------------------------------------------------------- + ir=0 + DO + r=ro+cosfac*radius1*(one-rbrack_l) + z=zo+sinfac*radius1*(one-rbrack_l) + CALL direct_get_bfield(r,z,bf,1) + dfdradius = -bf%psir*cosfac-bf%psiz*sinfac + dradius1 = -(psi0-bf%psi)/dfdradius + IF(dradius1<0.0)THEN + rbrack_l=rbrack_l*10 + IF(rbrack_l>radius1)THEN + CALL program_stop("find_fl_surface2 failed.") + ENDIF + ELSE + EXIT + ENDIF + + ir=ir+1 + IF (ir > direct_infinite_loop_count) THEN + direct_infinite_loop_flag = .TRUE. + CALL program_stop("Took too many steps to get rbrack_l.") + ENDIF + ENDDO +c----------------------------------------------------------------------- +c finding upper bound outside chosen surface +c dradius3 should always be less than 0. +c----------------------------------------------------------------------- + ir=0 + DO + r=ro+cosfac*radius1*(one+rbrack_h) + z=zo+sinfac*radius1*(one+rbrack_h) + CALL direct_get_bfield(r,z,bf,1) + dfdradius = -bf%psir*cosfac-bf%psiz*sinfac + dradius3 = -(psi0-bf%psi)/dfdradius + IF(dradius3>0.0)THEN + rbrack_h=rbrack_h*10 + IF(rbrack_h>radius1)THEN + CALL program_stop("find_fl_surface2 failed.") + ENDIF + ELSE + EXIT + ENDIF + + ir=ir+1 + IF (ir > direct_infinite_loop_count) THEN + direct_infinite_loop_flag = .TRUE. + CALL program_stop("Took too many steps to get rbrack_h.") + ENDIF + ENDDO +c----------------------------------------------------------------------- +c looping to narrow radii-intervals +c----------------------------------------------------------------------- + radii(1)=radius1-rbrack_l + radii(3)=radius1+rbrack_h + + ir=0 + DO + radii(2) = 0.5*(radii(3)+radii(1)) + IF(0.5*(radii(3)-radii(1)) < eps*ro)EXIT + + r=ro+cosfac*radii(1) + z=zo+sinfac*radii(1) + CALL direct_get_bfield(r,z,bf,1) + dfdradius = -bf%psir*cosfac-bf%psiz*sinfac + dradius1 = -(psi0-bf%psi)/dfdradius + + r=ro+cosfac*radii(2) + z=zo+sinfac*radii(2) + CALL direct_get_bfield(r,z,bf,1) + dfdradius = -bf%psir*cosfac-bf%psiz*sinfac + dradius2 = -(psi0-bf%psi)/dfdradius + + r=ro+cosfac*radii(3) + z=zo+sinfac*radii(3) + CALL direct_get_bfield(r,z,bf,1) + dfdradius = -bf%psir*cosfac-bf%psiz*sinfac + dradius3 = -(psi0-bf%psi)/dfdradius + + IF(dradius3>0.0 .OR. dradius1<0.0) + $ CALL program_stop("find_fl_surface2 error, needs debug.") + + IF(dradius2<0.0)THEN + radii(3)=radii(2) + ELSE + radii(1)=radii(2) + ENDIF + + ir=ir+1 + IF (ir > direct_infinite_loop_count) THEN + direct_infinite_loop_flag = .TRUE. + CALL program_stop("find_fl_surface2 failed.") + ENDIF + ENDDO +c----------------------------------------------------------------------- +c calculating r, z. +c----------------------------------------------------------------------- + r=ro+cosfac*radii(2) + z=zo+sinfac*radii(2) +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE find_fl_surface2 +c----------------------------------------------------------------------- +c subprogram 10. direct_xpt_y_out. +c makes y_out, splicing in analytic integrals. +c called when the numerical integral in direct_fl_int was terminated +c due to proximity to an x-point. can account for up to 2 x-points, +c if two are present, they must occur either side of the interval +c eta = [0.9pi,pi] for this script to work... aka they can't lie too +c close to the inboard midplane. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE direct_xpt_y_out(y_out1,len_y_out1,psifac,y_out, + $ len_y_out,eval_BpOnBt,eval_dy1,eval_xpt_tol,ix,diffs_raw, + $ outmat2f) + + REAL(r8), DIMENSION(0:,0:), INTENT(IN) :: y_out1 + REAL(r8),DIMENSION(1:4),INTENT(OUT) :: diffs_raw + REAL(r8),DIMENSION(1:17),INTENT(OUT) :: outmat2f + INTEGER, INTENT(IN) :: len_y_out1 + INTEGER, DIMENSION(2), INTENT(OUT) :: ix + REAL(r8), DIMENSION(0:,0:),INTENT(INOUT) :: + $ y_out + REAL(r8), INTENT(IN) :: psifac + REAL(r8), INTENT(OUT) :: eval_BpOnBt,eval_dy1,eval_xpt_tol + INTEGER,INTENT(INOUT) :: len_y_out + REAL(r8), DIMENSION(0:(4*nstepd+2*nstep2+6),0:4) :: y_outf + REAL(r8), DIMENSION(0:(nstepd+2),0:4) :: y_out2,y_out3, + $ y_out4 + REAL(r8), DIMENSION(0:(nstep2+2),0:4) :: yi1!,yi2 + TYPE(direct_bfield_type) :: bf + + INTEGER :: len_y_out2,len_y_out3,len_y_out4,len_y_outf + INTEGER :: ix1,ix2,i + ix=-1 +c----------------------------------------------------------------------- +c setting bounardies on integration interval: +c----------------------------------------------------------------------- + xpt_brackets(1,1)=y_out1(len_y_out1,0) +c----------------------------------------------------------------------- +c integrating opposite direction (clockwise). +c----------------------------------------------------------------------- + CALL direct_fl_int(psifac,twopi+0.1*pi,zero,y_out2, + $ len_y_out2,eval_BpOnBt,eval_dy1,eval_xpt_tol,.TRUE.,bf,.FALSE.) + IF(debug_xpt)WRITE(*,*)"Clockwise integral w. length:",len_y_out2 + IF(debug_xpt)WRITE(*,*)y_out2(len_y_out2-3:len_y_out2,0) + IF(debug_xpt)WRITE(*,*)y_out2(len_y_out2-3:len_y_out2,1) + IF(debug_xpt)WRITE(*,*)y_out2(len_y_out2-3:len_y_out2,2) + IF(debug_xpt)WRITE(*,*)y_out2(len_y_out2-3:len_y_out2,3) + IF(debug_xpt)WRITE(*,*)y_out2(len_y_out2-3:len_y_out2,4) +c----------------------------------------------------------------------- +c flagging case where one integral passed x-point while other +c did not, which should not happen. +c----------------------------------------------------------------------- + IF(y_out2(len_y_out2,0)y_out2(len_y_out2,0))THEN + CALL program_stop("Xpt int. err2, needs debug") + ENDIF + xpt_brackets(1,2)=y_out3(len_y_out3,0) + xpt_brackets(2,1)=y_out4(len_y_out4,0) +c----------------------------------------------------------------------- +c combining the two inboard numerical integrals into one +c anti-clockwise integral y_out +c----------------------------------------------------------------------- + IF(debug_xpt)WRITE(*,*)"y_out3 bounds:" + IF(debug_xpt)WRITE(*,*)y_out3(0,0),y_out3(len_y_out3,0) + IF(debug_xpt)WRITE(*,*)"y_out4 bounds:" + IF(debug_xpt)WRITE(*,*)y_out4(0,0),y_out4(len_y_out4,0) + CALL patch_two_a(y_out4,y_out3,len_y_out4,len_y_out3, + $ y_out,len_y_out) + IF(debug_xpt)WRITE(*,*)"y_out bounds:" + IF(debug_xpt)WRITE(*,*)y_out(0,0),y_out(len_y_out,0) + ix2=len_y_out +c----------------------------------------------------------------------- +c writing 1st analytic x-point integral and combining with inboard +c integrals. +c y_out, yi1 and y_outf all used as intermediaries in the splicing, +c but final result is in y_outf. +c----------------------------------------------------------------------- + CALL analytic_y_out(psifac,y_out4(len_y_out4,0), + $ y_out4(len_y_out4,2),y_out2(len_y_out2,0),y_out2(len_y_out2,2), + $ nstep2,yi1,diffs_raw,outmat2f) + IF(debug_xpt)WRITE(*,*)"yi1 bounds:" + IF(debug_xpt)WRITE(*,*)yi1(0,0),yi1(nstep2,0) + CALL patch_two_b(y_out,yi1,len_y_out,nstep2, + $ y_outf,len_y_outf) + + IF(debug_xpt)WRITE(*,*)"y_outf bounds:" + IF(debug_xpt)WRITE(*,*)y_outf(0,0),y_outf(len_y_outf,0) +c----------------------------------------------------------------------- +c writing 2nd analytic x-point integral and combining with y_outf. +c y_out, yi1 and y_outf all used as intermediaries in the splicing, +c but final result is in y_out. +c----------------------------------------------------------------------- + CALL analytic_y_out(psifac,y_out1(len_y_out1,0), + $ y_out1(len_y_out1,2),y_out3(len_y_out3,0),y_out3(len_y_out3,2), + $ nstep2,yi1,diffs_raw,outmat2f) + IF(debug_xpt)WRITE(*,*)"yi1 bounds: ROUND 2" + IF(debug_xpt)WRITE(*,*)yi1(0,0),yi1(nstep2,0) + CALL patch_two_b(yi1,y_outf,nstep2,len_y_outf, + $ y_out,len_y_out) + IF(debug_xpt)WRITE(*,*)"y_out bounds:" + IF(debug_xpt)WRITE(*,*)y_out(0,0),y_out(len_y_out,0) + ix2=ix2+nstep2 +c----------------------------------------------------------------------- +c reconstructing full 2pi integral, combining initial anti-clockwise +c numerical integral (y_out1) with y_out (inboard integrals plus +c both analytic x-point integrals), and then combining that with the +c second clockwise numerical integral (y_out2). +c----------------------------------------------------------------------- + IF(debug_xpt)WRITE(*,*)"y_out1 bounds:" + IF(debug_xpt)WRITE(*,*)y_out1(0,0),y_out1(len_y_out1,0) + CALL patch_two_b(y_out1,y_out,len_y_out1,len_y_out, + $ y_outf,len_y_outf) + IF(debug_xpt)WRITE(*,*)"y_outf bounds:" + IF(debug_xpt)WRITE(*,*)y_outf(0,0),y_outf(len_y_outf,0) + ix2=ix2+len_y_out+1 + CALL patch_two_c(y_outf,y_out2,len_y_outf,len_y_out2, + $ y_out,len_y_out,ix1) + ix2=ix2+nstep2+ix1+1 + ix(1)=ix1 + ix(2)=ix2 + ENDIF + IF(debug_xpt)WRITE(*,*)y_out(1:100,0) + IF(debug_xpt)WRITE(*,*)y_out(1:4,1) + IF(debug_xpt)WRITE(*,*)y_out(1:4,2) + IF(debug_xpt)WRITE(*,*)y_out(1:4,3) + IF(debug_xpt)WRITE(*,*)y_out(1:4,4) +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE direct_xpt_y_out +c----------------------------------------------------------------------- +c subprogram 11. direct_local_xpoint. +c finds location of nearby x-point where |Bp|=0 using Newton method. +c can and will search outside separatrix. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE direct_local_xpoint(r,z) + + REAL(r8), INTENT(INOUT) :: r,z + REAL(r8) :: ajac(2,2),det,dr,dz + REAL(r8), PARAMETER :: eps=1e-13 + TYPE(direct_bfield_type) :: bf + INTEGER :: ir +c----------------------------------------------------------------------- +c use newton iteration to find x-point. +c----------------------------------------------------------------------- + ir = 0 + DO + CALL direct_get_bfield(r,z,bf,2) + ajac(1,1)=bf%brr + ajac(1,2)=bf%brz + ajac(2,1)=bf%bzr + ajac(2,2)=bf%bzz + det=ajac(1,1)*ajac(2,2)-ajac(1,2)*ajac(2,1) + dr=(ajac(1,2)*bf%bz-ajac(2,2)*bf%br)/det + dz=(ajac(2,1)*bf%br-ajac(1,1)*bf%bz)/det + r=r+dr + z=z+dz + IF(ABS(dr) <= eps*r .AND. ABS(dz) <= eps*r)EXIT + + ir = ir+1 + IF (ir > direct_infinite_loop_count) THEN + direct_infinite_loop_flag = .TRUE. + CALL program_stop("Took too many steps to find x-point.") + ENDIF + ENDDO +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE direct_local_xpoint +c----------------------------------------------------------------------- +c subprogram 12. direct_saddle_angle. +c finds angle location of nearby saddle-node eigenvector using +c a binary search type algorithm. can find point where Bnu = 0 or +c Brho = 0. A newton method was attempted but generally failed to +c converge. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE direct_saddle_angle(rx,zx,rho,nustart_in,nu_var_in,nu, + $ Bcase,debug) + + REAL(r8), INTENT(IN) :: rx,zx,rho,nustart_in,nu_var_in + CHARACTER, INTENT(IN) :: Bcase + LOGICAL, INTENT(IN) :: debug + REAL(r8), INTENT(OUT) :: nu + + INTEGER, PARAMETER :: ird=4 + REAL(r8), PARAMETER :: nu_eps=1e-13 + + REAL(r8), DIMENSION(ird) :: nus,nu_tmp + REAL(r8) :: bnorm,Bout,nustep,pos,nustart,nu_var + INTEGER :: i,inuh + + + nustart=nustart_in + nu_var=nu_var_in + pos=one + inuh=0 +c----------------------------------------------------------------------- +c looping over narrower nu-intervals +c----------------------------------------------------------------------- + CALL direct_Blocal(rx,zx,nustart_in,rho,Bcase,bnorm) + IF(bnorm 80)THEN + direct_infinite_loop_flag = .TRUE. + CALL program_stop("Took too many steps to find x-pt angle.") + ENDIF + ENDDO +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE direct_saddle_angle +c----------------------------------------------------------------------- +c subprogram 13. direct_psisaddle. +c calculates the linear term of psi_in at the saddle point, as well +c as gamma, and vartheta (please refer to +c https://doi.org/10.1088/1361-6587/add9ca for definitions +c of these parameters). +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE direct_psisaddle(rx,zx, + $ nu,c11,gamma,vartheta,lincheck) + + REAL(r8), INTENT(IN) :: rx,zx + REAL(r8), DIMENSION(2), INTENT(IN) :: nu + REAL(r8), INTENT(OUT) :: c11,lincheck,gamma,vartheta + + REAL(r8), PARAMETER :: nuh_eps=1e-13, nuh_eps2=1e-6 + INTEGER :: ir + REAL(r8) :: nuh,nuperf,cosfac,sinfac,cosfact,sinfact + REAL(r8) :: psix,psinuh,r,Rlocal,Zlocal,x,y,chi + REAL(r8), DIMENSION(4) :: r_eps + TYPE(direct_bfield_type) :: bf + LOGICAL :: debug=.FALSE. + + r_eps(1) = 1e-7 + r_eps(2) = 1e-8 + r_eps(3) = 1e-9 + r_eps(4) = 1e-10 + ir=1 +c----------------------------------------------------------------------- +c finding correct initialisation point for nuh (read nu-half). +c gamma should be less than pi for a real x-point +c----------------------------------------------------------------------- + nuperf = (nu(1)+nu(2))/2 + gamma = nu(1)-nu(2) + nuh = nuperf + + IF (gamma > pi) THEN + CALL program_stop("Angle between separatrix legs > pi.") + ENDIF +c----------------------------------------------------------------------- +c finding angle where Brho = 0. +c----------------------------------------------------------------------- + CALL direct_saddle_angle(rx,zx,rx*r_eps(ir),nu(2),nu(1)-nu(2) + $ ,nuh,'r',.FALSE.) +c----------------------------------------------------------------------- +c make sure nuh is in [0,2pi). works even if nuh is negative. +c----------------------------------------------------------------------- + nuh = nuh - twopi*floor(nuh/twopi) +c----------------------------------------------------------------------- +c in the linear approaximation of psi at the x-point, Brho=0 halfway +c between the two separatrix legs. we check if we are within +c nuh_eps2 of this case +c----------------------------------------------------------------------- + lincheck = abs(nuh-nuperf) + IF(debug .OR. lincheck > nuh_eps2)THEN + PRINT"(A)","nu value where Brho=0 deviates from linear case by" + PRINT "(es16.10)", lincheck + ENDIF + IF(lincheck > nuh_eps2)THEN + CALL program_stop("psi isn't well approximated at the x-point") + ENDIF +c----------------------------------------------------------------------- +c defining rotated saddle-point coordinate frame to extract linear +c component. +c----------------------------------------------------------------------- + vartheta = nu(1)-pi/2.0 + vartheta = vartheta - twopi*floor(vartheta/twopi) + + r = r_eps(ir)*rx + cosfac=COS(nuh) + sinfac=SIN(nuh) + cosfact=COS(vartheta) + sinfact=SIN(vartheta) + + Rlocal = r*cosfac + Zlocal = r*sinfac + x = cosfact*Rlocal + sinfact*Zlocal + y = -sinfact*Rlocal + cosfact*Zlocal + chi = -COS(gamma)*x+SIN(gamma)*y +c----------------------------------------------------------------------- +c extracting linear component c11, where psi = psi(rx,zx)+c11*x*chi. +c----------------------------------------------------------------------- + CALL direct_get_bfield(rx,zx,bf,1) + psix = bf%psi + CALL direct_get_bfield(rx+r*cosfac,zx+r*sinfac,bf,1) + psinuh = bf%psi + + c11 = -(psinuh-psix)/(x*chi) +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE direct_psisaddle +c----------------------------------------------------------------------- +c subprogram 14. direct_xpoint. +c finds location and angles of nearby x-point, checks if it's inside +C separatrix. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE direct_xpoint(rin,zin,x_i,new_xpt) + + REAL(r8), INTENT(IN) :: rin,zin + INTEGER, INTENT(IN) :: x_i + LOGICAL, INTENT(OUT) :: new_xpt + + INTEGER, PARAMETER :: ird=4 + INTEGER :: i + REAL(r8), PARAMETER :: psi_eps=1e-4, r_eps1=1e-9 + REAL(r8) :: r,z + REAL(r8) :: c11,lincheck,gamma,vartheta + REAL(r8), DIMENSION(1:2) :: nu + REAL(r8) :: oangle,nu_var,Bnua,Bnub,Bnuc + TYPE(direct_bfield_type) :: bf + LOGICAL :: test_direct_local_xpoint,test_direct_saddle_angle + CHARACTER(64) :: message + + r=rin + z=zin + nu=0.0 + test_direct_local_xpoint=.FALSE. + test_direct_saddle_angle=.FALSE. + new_xpt=.FALSE. +c----------------------------------------------------------------------- +c testing direct_local_xpoint. +c----------------------------------------------------------------------- + IF(test_direct_local_xpoint)THEN + PRINT "(A)", "__________________________________________" + PRINT "(A)", "direct_local_xpoint output =>" + PRINT "(A)", "X-point number:" + PRINT "(i6)", x_i + PRINT "(A)", "eta:" + PRINT "(f16.10)", ATAN2(z-zo,r-ro) + CALL direct_get_bfield(ro,zo,bf,1) + PRINT "(A)", "psi at origin:" + PRINT "(f20.14)", bf%psi + PRINT "(A)", "r,z before direct_local_xpoint:" + PRINT "(f16.10)", r + PRINT "(f16.10)", z + CALL direct_get_bfield(r,z,bf,1) + PRINT "(A)", "|Bp| before direct_local_xpoint:" + PRINT "(f20.14)", SQRT(bf%br**2+bf%bz**2) + PRINT "(A)", "psi before direct_local_xpoint:" + PRINT "(f20.14)", bf%psi + CALL direct_local_xpoint(r,z) + PRINT "(A)", "r,z after direct_local_xpoint:" + PRINT "(f16.10)", r + PRINT "(f16.10)", z + CALL direct_get_bfield(r,z,bf,1) + PRINT "(A)", "|Bp| after direct_local_xpoint:" + PRINT "(f20.14)", SQRT(bf%br**2+bf%bz**2) + PRINT "(A)", "psi after direct_local_xpoint:" + PRINT "(f20.14)", bf%psi + !PRINT "(A)", "direct_local_xpoint output <=" + PRINT "(A)", "------------------------------------------" + ENDIF +c----------------------------------------------------------------------- +c finds x-point +c----------------------------------------------------------------------- + CALL direct_local_xpoint(r,z) + CALL direct_get_bfield(r,z,bf,1) +c----------------------------------------------------------------------- +c checks if this xpoint has already been identified +c----------------------------------------------------------------------- + IF(x_i>0)THEN + DO i=1,(x_i-1),+1 + IF(ABS(r-rxs(i))<(0.0001*ro) .AND. + $ ABS(z-zxs(i))<(0.0001*ro))THEN + RETURN + ENDIF + ENDDO + ENDIF +c----------------------------------------------------------------------- +c fills out global module variables +c----------------------------------------------------------------------- + rxs(x_i) = r + zxs(x_i) = z +c----------------------------------------------------------------------- +c checks if outside separatrix. +c----------------------------------------------------------------------- + IF (bf%psi < -psi_eps*psio) THEN + outside_sep(x_i)=.TRUE. + ELSE + outside_sep(x_i)=.FALSE. + ENDIF +c----------------------------------------------------------------------- +c updating xpt_etas with a more exact value. +c----------------------------------------------------------------------- + xpt_etas(x_i)=ATAN2(zxs(x_i)-zo,rxs(x_i)-ro) + xpt_etas(x_i)=xpt_etas(x_i) - twopi*floor(xpt_etas(x_i)/twopi) +c----------------------------------------------------------------------- +c finding angle of o-point from perspective of x-point (oangle). +c----------------------------------------------------------------------- + oangle = xpt_etas(x_i)+pi + oangle = oangle - twopi*floor(oangle/twopi) + nu_var = pi/2 +c----------------------------------------------------------------------- +c calculating separatrix leg angles (nu). nu(1) is always more than +c oangle, nu(2) is always less than oangle. +c----------------------------------------------------------------------- + CALL direct_saddle_angle(rxs(x_i),zxs(x_i),r_eps1*rxs(x_i),oangle, + $ nu_var,nu(1),'n',.FALSE.) + CALL direct_saddle_angle(rxs(x_i),zxs(x_i),r_eps1*rxs(x_i),oangle, + $ -nu_var,nu(2),'n',.FALSE.) +c----------------------------------------------------------------------- +c testing direct_saddle_angle +c----------------------------------------------------------------------- + IF(test_direct_saddle_angle)THEN + CALL direct_Blocal(rxs(x_i),zxs(x_i),nu(1),r_eps1*rxs(x_i), + $ 'n',Bnua) + CALL direct_Blocal(rxs(x_i),zxs(x_i),oangle,r_eps1*rxs(x_i), + $ 'n',Bnub) + CALL direct_Blocal(rxs(x_i),zxs(x_i),nu(2),r_eps1*rxs(x_i), + $ 'n',Bnuc) + PRINT "(A)", "__________________________________________" + PRINT "(A)", "direct_saddle_angle output =>" + PRINT "(A)", "first x-point leg's angle nu:" + PRINT "(f17.14)", nu(1)/pi + PRINT "(A)", "first leg B_nu:" + PRINT "(es16.3)", Bnua + PRINT "(A)", "nu angle pointing from x-pt to mag. axis:" + PRINT "(f17.14)", oangle/pi + PRINT "(A)", "B_nu at angle pointing from x-pt to mag. axis:" + PRINT "(es16.3)", Bnub + PRINT "(A)", "second x-point leg's angle nu:" + PRINT "(f17.14)", nu(2)/pi + PRINT "(A)", "second leg B_nu:" + PRINT "(es16.3)", Bnuc + PRINT "(A)", "angle between x-point legs (gamma):" + PRINT "(f17.14)", (nu(1)-nu(2))/pi + PRINT "(A)", "------------------------------------------" + ENDIF +c----------------------------------------------------------------------- +c calculating xpoint angles gamma, vartheta, and linear term c11. +c----------------------------------------------------------------------- + CALL direct_psisaddle(r,z,nu,c11,gamma,vartheta,lincheck) +c----------------------------------------------------------------------- +c filling out x-point module variables. +c----------------------------------------------------------------------- + xpt_c11s(x_i) = c11 + xpt_gammas(x_i) = gamma + xpt_varthetas(x_i) = vartheta +c----------------------------------------------------------------------- +c acknowledges new xpoint after all module variables filled +c----------------------------------------------------------------------- + new_xpt=.TRUE. +c----------------------------------------------------------------------- +c regular print statements. eta bracket describes the region where +c the analytic formulas will take over from the numerical +c integrator. +c----------------------------------------------------------------------- +c501 FORMAT(1x,i1," x-point(s) detected:") +c502 FORMAT(1x,"over eta bracket (",f10.8,",",f10.8,")") +600 FORMAT(1x," x-point detected at eta, r = ",f6.4,",",f6.4) +601 FORMAT(1x," 2nd x-point detected at eta, r = ",f6.4,",",f6.4) +603 FORMAT(1x,"integrals applied for eta = (", + $ f7.5,",",f7.5,").") + IF(verbose)THEN + PRINT "(A)", " _______________________________________________" + IF(x_i==1)THEN + WRITE(message,600)xpt_etas(x_i), + $ SQRT((zxs(x_i)-zo)**2+(rxs(x_i)-ro)**2) + PRINT "(A)",message + ELSEIF(x_i==2)THEN + WRITE(message,601)xpt_etas(x_i), + $ SQRT((zxs(x_i)-zo)**2+(rxs(x_i)-ro)**2) + PRINT "(A)",message + ENDIF + !WRITE(message,603)xpt_brackets(x_i,1),xpt_brackets(x_i,2) + !PRINT "(A)",message + PRINT "(A)", " -----------------------------------------------" + ENDIF +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE direct_xpoint +c----------------------------------------------------------------------- +c subprogram 15. direct_saddle_coords. +c transforms us into a local coordinate frame aligned with one leg +c of the x-point. takes x-point data from module variables. +c r1 is radius from the magnetic axis, not major radius +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE direct_saddle_coords(x_i,r1,eta1,x,y,chi,usevth2) + + INTEGER, INTENT(IN) :: x_i + LOGICAL, INTENT(IN) :: usevth2 + REAL(r8), INTENT(IN) :: r1,eta1 + REAL(r8), INTENT(OUT) :: x,y,chi + + REAL(r8) :: cosfac,sinfac,cosfact,sinfact + REAL(r8) :: Rshft,Zshft,Rlocal,Zlocal + +c----------------------------------------------------------------------- +c precalculating trig components. +c----------------------------------------------------------------------- + cosfac=COS(eta1) + sinfac=SIN(eta1) + IF(usevth2)THEN + cosfact=COS(xpt_varthetas2(x_i)) + sinfact=SIN(xpt_varthetas2(x_i)) + ELSE + cosfact=COS(xpt_varthetas(x_i)) + sinfact=SIN(xpt_varthetas(x_i)) + ENDIF +c----------------------------------------------------------------------- +c local R,Z coordinates centered on x-point. +c----------------------------------------------------------------------- + Rshft = ro + r1*cosfac + Zshft = zo + r1*sinfac + + Rlocal = Rshft-rxs(x_i) + Zlocal = Zshft-zxs(x_i) +c----------------------------------------------------------------------- +c rotating local R,Z coordinates to align with x-point leg that +c approaches the x-point by travelling anticlockwise around the +c separatrix. this is captured in the vartheta variable +c----------------------------------------------------------------------- + x = cosfact*Rlocal + sinfact*Zlocal + y = -sinfact*Rlocal + cosfact*Zlocal +c----------------------------------------------------------------------- +c calculating chi angle variable, defined such that nabla chi is +c orthogonal to the x-point leg that leaves the +c x-point when travelling anticlockwise around the separatrix. +c----------------------------------------------------------------------- + IF(usevth2)THEN + chi = -COS(xpt_gammas2(x_i))*x+SIN(xpt_gammas2(x_i))*y + ELSE + chi = -COS(xpt_gammas(x_i))*x+SIN(xpt_gammas(x_i))*y + ENDIF +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE direct_saddle_coords +c----------------------------------------------------------------------- +c subprogram 16. direct_saddle_coords_inv. +c inverse of direct_saddle_coords. takes x,y, returns R,Z,rho,eta. +c rho is minor radius. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE direct_saddle_coords_inv(x_i,x,y,R,Z,rho,eta,usevth2) + + INTEGER, INTENT(IN) :: x_i + LOGICAL, INTENT(IN) :: usevth2 + REAL(r8), INTENT(IN) :: x,y + REAL(r8), INTENT(OUT) :: R,Z,rho,eta + + REAL(r8) :: cosfact,sinfact + REAL(r8) :: Rlocal,Zlocal +c----------------------------------------------------------------------- +c precalculating trig components. +c----------------------------------------------------------------------- + IF(usevth2)THEN + cosfact=COS(xpt_varthetas2(x_i)) + sinfact=SIN(xpt_varthetas2(x_i)) + ELSE + cosfact=COS(xpt_varthetas(x_i)) + sinfact=SIN(xpt_varthetas(x_i)) + ENDIF +c----------------------------------------------------------------------- +c inverse coordinate transform. +c----------------------------------------------------------------------- + Rlocal = cosfact*x - sinfact*y + Zlocal = sinfact*x + cosfact*y + + R = Rlocal + rxs(x_i) + Z = Zlocal + zxs(x_i) + + rho = SQRT((R-ro)**2+(Z-zo)**2) + eta = ATAN2(Z-zo,R-ro) +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE direct_saddle_coords_inv +c----------------------------------------------------------------------- +c subprogram 17. int_validity_check. +c checks if we are close enough to the x-point for the analytic +c formulas to be valid. This is Eq. 12 in +c https://doi.org/10.1088/1361-6587/add9ca. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE int_validity_check(r,z,r_loc,eta,valid,eval_xpt_tol) + + REAL(r8), INTENT(IN) :: r,z,r_loc + REAL(r8), INTENT(INOUT) :: eta + REAL(r8), DIMENSION(1:2) :: xpt_dists + LOGICAL, INTENT(OUT) :: valid + REAL(r8), INTENT(OUT) :: eval_xpt_tol + REAL(r8) :: x,y,chi + INTEGER :: ixpt + !x1_tol gonna +c----------------------------------------------------------------------- +c find nearest x-point. +c----------------------------------------------------------------------- + IF(num_xpts==0)CALL program_stop("X-point init missed something.") + xpt_dists(1)=SQRT((r-rxs(1))**2+(z-zxs(1))**2) + IF(num_xpts==2)THEN + xpt_dists(2)=SQRT((r-rxs(2))**2+(z-zxs(2))**2) + ELSE + xpt_dists(2)=1d99 + ENDIF + ixpt=MINLOC(xpt_dists,1) +c----------------------------------------------------------------------- +c check we are actually near this x-point in angle-space. +c----------------------------------------------------------------------- + eta=eta-twopi*floor(eta/twopi) + IF(ABS(xpt_etas(ixpt)-eta)>0.1*pi.OR.xpt_dists(ixpt)>0.1*ro)THEN + CALL program_stop("Poloidal field vanishing far from x-point") + ENDIF +c----------------------------------------------------------------------- +c if etaABS(chi))CALL program_stop("angle sanity ch. fail 1") + !^We are approaching anticlockwise, travelling down the x-axis. + !The int. point should be closer to the first leg of the sep. + !then the second + eval_xpt_tol=ABS(x/xpt_dists(ixpt)) + ELSE + IF(ABS(chi)>ABS(x))CALL program_stop("angle sanity ch. fail 2") + eval_xpt_tol=ABS(chi/xpt_dists(ixpt)) + !chi is the orthogonal distance from the second leg of the + !x-point. The second leg is the one pointing in the + !anticlockwise direction as we travel and the psi surface in + !a normal poloidal cross section image. + ENDIF + valid=eval_xpt_tol", valid +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE int_validity_check +c----------------------------------------------------------------------- +c subprogram 18. direct_analytic_ints. +c takes the linear term c11, angle and position data, and calculates +c integrals y(1), y(2), y(3), y(4) replacing direct_fl_int. assumes +c logarithmic trajectories near the saddle point. +c for more information see https://doi.org/10.1088/1361-6587/add9ca. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE direct_analytic_ints(x_i,r1,eta1,eta2,yi1,yi2,yi3,yi4, + $ outmat,outmat2,debug) + + REAL(r8), INTENT(IN) :: r1,eta1,eta2 + INTEGER, INTENT(IN) :: x_i + LOGICAL, INTENT(IN) :: debug + LOGICAL :: usevth2 + + REAL(r8), INTENT(OUT) :: yi1,yi2,yi3,yi4 + REAL(r8), DIMENSION(1:4,1:2,0:1),INTENT(OUT) :: outmat + REAL(r8), DIMENSION(1:8),INTENT(OUT) :: outmat2 + + REAL(r8) :: x1,y1,x2,y2,chi,J1,J2,J3,bt,b,root,xo,yo,chio,d,C00 + REAL(r8) :: y3d,y3r1,y3r2,singmt,a,c,yi2b,yi3b,cosgam + REAL(r8) :: cosvt,sinvt,cotgam,singam,cscgam,etax,vartheta,gamma + TYPE(direct_bfield_type) :: bf + + outmat=0.0 + usevth2=.TRUE. +c----------------------------------------------------------------------- +c getting local coordinates. +c----------------------------------------------------------------------- + CALL direct_saddle_coords(x_i,r1,eta1,x1,y1,chi,.TRUE.) + CALL direct_saddle_coords(x_i,zero,zero,xo,yo,chio,.TRUE.) +c----------------------------------------------------------------------- +c pre-calculationg useful terms +c----------------------------------------------------------------------- + IF(usevth2)THEN + vartheta=xpt_varthetas2(x_i) + gamma=xpt_gammas2(x_i) + ELSE + vartheta=xpt_varthetas(x_i) + gamma=xpt_gammas(x_i) + ENDIF + + cosvt = COS(vartheta) + sinvt = SIN(vartheta) + + cotgam = one/TAN(gamma) + singam = SIN(gamma) + cosgam = COS(gamma) + singmt = SIN(gamma-vartheta) + + cscgam = one/singam + etax = xpt_etas(x_i) + CALL direct_get_bfield(rxs(x_i),zxs(x_i),bf,1) + bt=abs(bf%f/rxs(x_i)) + b=bt +c----------------------------------------------------------------------- +c calculating x2 by assuming logarithmic Bp trajectory near x-point. +c----------------------------------------------------------------------- + root = -sqrt((xo*TAN(eta2-vartheta)-yo)**2 + $-4*one*(cotgam-TAN(eta2-vartheta))*(x1*y1-x1*x1*cotgam)) + x2 = (one/2)*(yo-xo*TAN(eta2-vartheta)+root)/ + $ (cotgam-TAN(eta2-vartheta)) + y2 = (x2-xo)*TAN(eta2-vartheta)+yo +c----------------------------------------------------------------------- +c outputting coordinate information. +c outmat(1,1,0) = x1, outmat(1,2,0) = y1 +c outmat(2,1,0) = R1, outmat(2,2,0) = Z1 +c outmat(3,1,0) = rho1, outmat(3,2,0) = eta1 +c outmat(1,1,1) = x2, outmat(1,2,1) = y2 +c outmat(2,1,1) = R2, outmat(2,2,1) = Z2 +c outmat(3,1,1) = rho2, outmat(3,2,1) = eta2 +c 1 and 2 refer to initial and final point respectively. rho is +c minor radius, eta poloidal angle. output yi2=rho2=outmat(3,1,1) +c----------------------------------------------------------------------- + IF(debug)THEN + outmat(1,1,0)=x1 + outmat(1,1,1)=x2 + outmat(1,2,0)=y1 + outmat(1,2,1)=y2 + + CALL direct_saddle_coords_inv(x_i,x1,y1,outmat(2,1,0), + $ outmat(2,2,0),outmat(3,1,0),outmat(3,2,0),.TRUE.) + ENDIF + CALL direct_saddle_coords_inv(x_i,x2,y2,outmat(2,1,1), + $ outmat(2,2,1),outmat(3,1,1),outmat(3,2,1),.TRUE.) +c----------------------------------------------------------------------- +c evaluating y1 +c----------------------------------------------------------------------- + yi1 = -(cscgam/xpt_c11s(x_i))*(rxs(x_i)*LOG(x2/x1) + $ +(x2-x1)*(cosvt-(one-x1/x2)*cotgam*sinvt) + $ -(one-x1/x2)*y1*sinvt) +c----------------------------------------------------------------------- +c evaluating y2 +c----------------------------------------------------------------------- + yi2=outmat(3,1,1) + + !The following formula for minor radius integral is wrong... + IF(debug)yi2b = -(cscgam/xpt_c11s(x_i))*abs(-SIN(etax+gamma + $ -vartheta)*(x2-x1) + $ +singam*SIN(etax-vartheta)*(y1-x1)*(one-x1/x2)) +c----------------------------------------------------------------------- +c evaluating y3 +c----------------------------------------------------------------------- + a=cscgam*singmt + c=x1*(x1*cotgam-y1)*sinvt + + y3d=rxs(x_i)**2.d0-4.d0*a*c + + y3r1=(one/(2*a))*(-rxs(x_i)+SQRT(y3d)) + y3r2=(one/(2*a))*(-rxs(x_i)-SQRT(y3d)) + + yi3 = -(cscgam/xpt_c11s(x_i))*(one/SQRT(y3d))* + $ LOG((x2-y3r1)*(x1-y3r2)/((x2-y3r2)*(x1-y3r1))) + IF(debug)yi3b = -(cscgam/(xpt_c11s(x_i)*rxs(x_i)))*(LOG(x2/x1) + $ -(cosvt-(one-x1/x2)*cotgam*sinvt)*(x2-x1)/rxs(x_i) + $ +sinvt*y1*(one-x1/x2)/rxs(x_i)) + + IF(debug)THEN + outmat(4,1,0)=yi2 + outmat(4,1,1)=yi2b + outmat(4,2,0)=yi3 + outmat(4,2,1)=yi3b + ENDIF +c----------------------------------------------------------------------- +c y4 special case. +c----------------------------------------------------------------------- + IF(abs(power_bp)>1e-13)THEN + IF(.NOT. debug)THEN + CALL program_stop("Analytic int incompat. w power_bp =/= 0") + ELSE + WRITE(*,*)"Analytic integrals incompatible w power_bp != 0" + WRITE(*,*)"Printing data" + ENDIF + !Refer to https://doi.org/10.1088/1361-6587/add9ca for the + !the full writeup of the y4 formula. + !Need all ingredients for these Appellf1 functions: + !Ingredients for J1, J2, J3 in Mathematica: + !d = x1*y1-x1*x1 + C00=-cscgam*(b**power_b)*(ABS(xpt_c11s(x_i)**(power_bp-one)))* + $ (xpt_c11s(x_i)/ABS(xpt_c11s(x_i))) !keeping sign consistent + outmat2(1)=y1 + outmat2(2)=x1 + outmat2(3)=x2 + outmat2(4)=gamma + outmat2(5)=power_bp + + outmat2(6)=C00* + $ rxs(x_i)**(one-power_r-power_bp) + outmat2(7)=C00* + $ (one-power_r-power_bp)* + $ (cosvt-cotgam*sinvt)/(rxs(x_i)**(power_r+power_bp)) + outmat2(8)=C00* + $ (-one)*(one-power_r-power_bp)* + $ sinvt*(x1*y1-x1*x1*cotgam)/(rxs(x_i)**(power_r+power_bp)) + + !yi4=outmat2(6)*J1+outmat2(7)*J2+outmat2(8)*J3 + yi4=0.0 + + !Full writeup: + !yi4 = -cscgam*(b**power_b)*(xpt_c11s(x_i)**(power_bp-one))*( + !$ J1*rxs(x_i)**(one-power_r-power_bp)+ + !$ J2*(one-power_r-power_bp)* + !$ (cosvt-cotgam*sinvt)/(rxs(x_i)**(power_r+power_bp)) + !$ -J3*(one-power_r-power_bp)* + !$ sinvt*(x1*y1-x1*x1)/(rxs(x_i)**(power_r+power_bp))) + + !AppellF1 entries (Can't appel package since xf1,yf1 complex) + !a(1)=-alpha_p/2.d0 + !b1(1)=-alpha_p/2.d0 + !b2(1)=-alpha_p/2.d0 + !c(1)=1.d0-alpha_p/2.d0 + + !a(2)=1.d0/2.d0-alpha_p/2.d0 + !b1(2)=-alpha_p/2.d0 + !b2(2)=-alpha_p/2.d0 + !c(2)=3.d0/2.d0-alpha_p/2.d0 + + !a(3)=-1.d0/2.d0-alpha_p/2.d0 + !b1(3)=-alpha_p/2.d0 + !b2(3)=-alpha_p/2.d0 + !c(3)=1.d0/2.d0-alpha_p/2.d0 + + !xf1=x**2/(d*cosgam*singam - sqrt(-d**2*singam**4)) + !yf1=xf1 + ELSE + outmat2(1)=0.d0 + outmat2(2)=0.d0 + outmat2(3)=0.d0 + outmat2(4)=0.d0 + outmat2(5)=0.d0 + outmat2(6)=0.d0 + outmat2(7)=0.d0 + outmat2(8)=0.d0 + + J1 = LOG(x2/x1) + J2 = x2-x1 + J3 = -(one/x2)+one/x1 + yi4 = -cscgam*(b**power_b)*(xpt_c11s(x_i)**(power_bp-one))*( + $ J1*rxs(x_i)**(one-power_r)+ + $ J2*(one-power_r)*(cosvt-cotgam*sinvt)/(rxs(x_i)**power_r) + $ -J3*(one-power_r)*sinvt*(x1*y1-x1*x1*cotgam)/(rxs(x_i)**power_r)) + ENDIF +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE direct_analytic_ints +c----------------------------------------------------------------------- +c subprogram 19. interp_through_zero. +c adds a spline knot at zero (or twopi). +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE interp_through_zero(y_out,len_y_out,i_zero) + + REAL(r8), DIMENSION(0:,0:), INTENT(INOUT) :: y_out + REAL(r8), DIMENSION(0:(4*nstepd+2*nstep2+6),0:4) :: y_cpy + INTEGER, INTENT(INOUT) :: len_y_out + INTEGER, INTENT(OUT) :: i_zero + REAL(r8), DIMENSION(0:4) :: yinterim1 + INTEGER :: i,len_y_out_tmp + LOGICAL :: interpd + !force_twopi=.TRUE. +c----------------------------------------------------------------------- +c if integrals pass through twopi, we add a point there. +c----------------------------------------------------------------------- + interpd=.FALSE. + + y_cpy(0:len_y_out,0)=y_out(0:len_y_out,0) + y_cpy(0:len_y_out,1)=y_out(0:len_y_out,1) + y_cpy(0:len_y_out,2)=y_out(0:len_y_out,2) + y_cpy(0:len_y_out,3)=y_out(0:len_y_out,3) + y_cpy(0:len_y_out,4)=y_out(0:len_y_out,4) + + len_y_out_tmp=len_y_out + DO i=0,len_y_out_tmp-1,+1 + IF(y_out(i+1,0)==twopi)THEN + WRITE(*,*)"patch_two_c no interp at 0 required, s2" + i_zero=i+1 + EXIT !no interp requied + ELSEIF(y_out(i,0)>twopi .AND. + $ y_out(i+1,0)etaiii>etaii>etai. Also assumes we aren't going more +c than twopi around the loop. Stitches splines together at +c 0.5*(etaii+etaiii). Returns y_out, len_y_out, which will seem +c like a single integral. +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE patch_two_a(y_out1,y_out2, + $ len_y_out1,len_y_out2,y_out,len_y_out) + + REAL(r8), DIMENSION(0:,0:), INTENT(INOUT) :: y_out1,y_out2 + REAL(r8), DIMENSION(0:,0:), INTENT(INOUT) :: y_out + REAL(r8), DIMENSION(0:(4*nstepd+2*nstep2+6),0:4) :: y_cpy + REAL(r8), DIMENSION(0:4) :: yinterim1,yinterim2 + INTEGER, INTENT(IN) :: len_y_out1,len_y_out2 + INTEGER, INTENT(INOUT) :: len_y_out + REAL(r8) :: etai,etaii,etaiii,etaiv,eta_s + INTEGER :: ii,i + y_out=0.0 + len_y_out=0 +c----------------------------------------------------------------------- +c assert angles. +c----------------------------------------------------------------------- + etaii =y_out1(0,0) + etaiv =y_out1(len_y_out1,0) + etaiii=y_out2(0,0) + etai =y_out2(len_y_out2,0) + IF(.NOT.((etaiv>etaiii).AND.(etaiii>etaii).AND.(etaii>etai)))THEN + WRITE(*,*)"Incorrect use of patch_two_a (s1)." + WRITE(*,*)y_out2(len_y_out2,0) + WRITE(*,*)y_out2(len_y_out2-1,0) + WRITE(*,*)y_out2(0,0) + WRITE(*,*)"Angles:" + WRITE(*,*)etaiv + WRITE(*,*)etaiii + WRITE(*,*)etaii + WRITE(*,*)etai + CALL program_stop("Incorrect use of patch_two_a. (s1)") + ELSEIF((etaiv-etai)>twopi)THEN + CALL program_stop("Incorrect use of patch_two_a. (s2)") + ENDIF +c----------------------------------------------------------------------- +c find main stitch point. +c----------------------------------------------------------------------- + eta_s=0.5*(etaiii+etaii) +c----------------------------------------------------------------------- +c create stich points. +c----------------------------------------------------------------------- + CALL yinterim(y_out1,len_y_out1,eta_s,yinterim1) + CALL yinterim(y_out2,len_y_out2,eta_s,yinterim2) +c----------------------------------------------------------------------- +c checking radii are similar. +c----------------------------------------------------------------------- + IF(ABS(yinterim1(2)-yinterim2(2))/ro > r_tol)THEN + WRITE(*,*)"patch_two_a:" + WRITE(*,*)ABS(yinterim1(2)-yinterim2(2))/ro + CALL program_stop( + $ "Radial discontinuity > r_tol during field line stitching.") + ENDIF +c----------------------------------------------------------------------- +c renormalising y_out1 so the integrals are zero at stitch point. +c----------------------------------------------------------------------- + DO i=0,len_y_out1,+1 + y_out1(i,1)=y_out1(i,1)-yinterim1(1) + y_out1(i,3)=y_out1(i,3)-yinterim1(3) + y_out1(i,4)=y_out1(i,4)-yinterim1(4) + ENDDO +c----------------------------------------------------------------------- +c resigning those integrals of absolute values in y_out2 +c----------------------------------------------------------------------- + y_cpy(0:len_y_out2,0)=y_out2(0:len_y_out2,0) + y_cpy(0:len_y_out2,1)=y_out2(0:len_y_out2,1) + y_cpy(0:len_y_out2,2)=y_out2(0:len_y_out2,2) + y_cpy(0:len_y_out2,3)=y_out2(0:len_y_out2,3) + y_cpy(0:len_y_out2,4)=y_out2(0:len_y_out2,4) + DO i=len_y_out2,0,-1 + y_out2(i,1)=-(y_cpy(i,1)-y_cpy(len_y_out2,1)) + y_out2(i,3)=-(y_cpy(i,3)-y_cpy(len_y_out2,3)) + y_out2(i,4)=-(y_cpy(i,4)-y_cpy(len_y_out2,4)) + ENDDO +c----------------------------------------------------------------------- +c writing y_out2 to y_out. +c----------------------------------------------------------------------- + ii=0 + DO i=len_y_out2,0,-1 + IF(y_out2(i,0)eta_s)THEN + !Theta and radius values just flat equals + y_out(ii,0)=y_out1(i,0) + y_out(ii,2)=y_out1(i,2) + + !Integrals need to add previous: + y_out(ii,1)=y_out1(i,1)+yinterim2(1) + y_out(ii,3)=y_out1(i,3)+yinterim2(3) + y_out(ii,4)=y_out1(i,4)+yinterim2(4) + ii=ii+1 + ENDIF + ENDDO + len_y_out=ii-1 +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE patch_two_a +c----------------------------------------------------------------------- +c subprogram 21. patch_two_b. +c combines two y_outs together. assumes y_out1 and y_out2 are both +c travelling anticlockwise and start/end at the same point. +c Returns y_out, len_y_out, which will seem like a single integral. +c Assumes integral doens't travel more than twopi in eta. +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE patch_two_b(y_out1,y_out2, + $ len_y_out1,len_y_out2,y_out,len_y_out) + + REAL(r8), DIMENSION(0:,0:), INTENT(IN) :: y_out1,y_out2 + REAL(r8), DIMENSION(0:,0:), INTENT(INOUT) :: y_out + INTEGER, INTENT(IN) :: len_y_out1,len_y_out2 + INTEGER, INTENT(INOUT) :: len_y_out + REAL(r8) :: etai,etaii,etaiii,etaiv + INTEGER :: ii,i + y_out=0.0 + len_y_out=0 +c----------------------------------------------------------------------- +c assert angles. +c----------------------------------------------------------------------- + etai =y_out1(0,0) + etaii =y_out1(len_y_out1,0) + etaiii=y_out2(0,0) + etaiv =y_out2(len_y_out2,0) + !IF(.NOT.(((etaiv>etaiii).AND.(etaiii==etaii).AND.(etaii>etai)).OR. + !$((etaiv>etaiii).AND.(etaiii>etaii).AND.((etaiii-etaii)<0.0001*pi) + !$.AND. (etaii>etai))))THEN + IF(.NOT.((etaiv>etaiii).AND.(etaiii==etaii).AND.(etaii>etai)))THEN + WRITE(*,*)"Incorrect use of patch_two_b (s1)." + WRITE(*,*)y_out2(len_y_out2,0) + WRITE(*,*)y_out2(len_y_out2-1,0) + WRITE(*,*)y_out2(0,0) + WRITE(*,*)"Angles:" + WRITE(*,*)etaiv + WRITE(*,*)etaiii + WRITE(*,*)etaii + WRITE(*,*)etai + CALL program_stop("Incorrect use of patch_two_b (s1).") + ELSEIF((etaiv-etai)>twopi)THEN + CALL program_stop("Incorrect use of patch_two_b (s2).") + ENDIF +c----------------------------------------------------------------------- +c checking radii are similar. +c----------------------------------------------------------------------- + IF(ABS(y_out1(len_y_out1,2)-y_out2(0,2))/ro > r_tol)THEN + WRITE(*,*)"patch_two_b:" + WRITE(*,*)y_out1(0,2) + WRITE(*,*)y_out2(0,2) + WRITE(*,*)ABS(y_out1(len_y_out1,2)-y_out2(0,2))/ro + CALL program_stop( + $ "Radial discontinuity > r_tol during field line stitching.") + ENDIF +c----------------------------------------------------------------------- +c adding y_out1 to y_out. +c----------------------------------------------------------------------- + ii=0 + DO i=0,len_y_out1,+1 + y_out(ii,:)=y_out1(i,:) + ii=ii+1 + ENDDO +c----------------------------------------------------------------------- +c writing y_out2 to y_out. +c----------------------------------------------------------------------- + DO i=1,len_y_out2,+1 + y_out(ii,0)=y_out2(i,0) + y_out(ii,2)=y_out2(i,2) + + y_out(ii,1)=y_out2(i,1)+y_out1(len_y_out1,1) + y_out(ii,3)=y_out2(i,3)+y_out1(len_y_out1,3) + y_out(ii,4)=y_out2(i,4)+y_out1(len_y_out1,4) + ii=ii+1 + ENDDO + len_y_out=ii-1 +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE patch_two_b +c----------------------------------------------------------------------- +c subprogram 22. patch_two_c. +c combines two y_outs together. assumes y_out1 is anticlockwise, +c from etai>=0 to etaii, and y_out2 is clockwise from +c etaiv>twopi to etaiii=etaii. +c iswtch is the index where the output spline y_out switches from +c y_out1 to y_out2 (technically index of final y_out1 point) +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE patch_two_c(y_out1,y_out2, + $ len_y_out1,len_y_out2,y_out,len_y_out,iswtch) + + REAL(r8), DIMENSION(0:,0:), INTENT(INOUT) :: y_out1,y_out2 + REAL(r8), DIMENSION(0:(4*nstepd+2*nstep2+6),0:4) :: y_cpy + REAL(r8), DIMENSION(0:,0:), INTENT(INOUT) :: y_out + REAL(r8), DIMENSION(0:4) :: yinterim1,yinterim2 + REAL(r8) :: etai,etaii,etaiii,etaiv,eta_s + INTEGER, INTENT(INOUT) :: len_y_out1,len_y_out2 + INTEGER, INTENT(INOUT) :: len_y_out + INTEGER, INTENT(OUT) :: iswtch + INTEGER :: ii,i1,i_zero,i + y_out=0.0 + len_y_out=0 +c----------------------------------------------------------------------- +c assert angles. +c----------------------------------------------------------------------- + etai =y_out1(0,0) + etaii =y_out1(len_y_out1,0) + etaiv =y_out2(0,0) + etaiii=y_out2(len_y_out2,0) + IF(.NOT.((etaiv>etaiii).AND.(etaiii==etaii).AND.(etaii>etai).AND. + $ (etai==zero).AND.(etaiv>twopi).AND.(etaiii r_tol)THEN + WRITE(*,*)"patch_two_c:" + WRITE(*,*)ABS(yinterim1(2)-yinterim2(2))/ro + CALL program_stop( + $ "Radial discontinuity > r_tol during field line stitching.") + ENDIF +c----------------------------------------------------------------------- +c renormalising y_out1 so the integrals are zero at stitch point. +c----------------------------------------------------------------------- + DO i=0,len_y_out1,+1 + y_out1(i,1)=y_out1(i,1)-yinterim1(1) + y_out1(i,3)=y_out1(i,3)-yinterim1(3) + y_out1(i,4)=y_out1(i,4)-yinterim1(4) + ENDDO +c----------------------------------------------------------------------- +c adding start of y_out2 (eta twopi to stitch point) to y_out. +c----------------------------------------------------------------------- + ii=0 + DO i=len_y_out2,0,-1 + IF(y_out2(i,0)>=twopi .AND. y_out2(i,0)(eta_s-twopi))THEN + y_out(ii,0)=y_out1(i,0) + y_out(ii,2)=y_out1(i,2) + + y_out(ii,1)=y_out1(i,1)+y_out(i1,1) + y_out(ii,3)=y_out1(i,3)+y_out(i1,3) + y_out(ii,4)=y_out1(i,4)+y_out(i1,4) + + ii=ii+1 + iswtch=iswtch+1 + ENDIF + ENDDO + i1=ii-1 + iswtch=i1 !end of y_out1 in y_out +c----------------------------------------------------------------------- +c writing rest of y_out2s to y_out. +c----------------------------------------------------------------------- + iswtch=0 + DO i=len_y_out2,0,-1 + IF(y_out2(i,0)<=twopi .AND. y_out2(i,0)>=y_out(i1,0))THEN + y_out(ii,0)=y_out2(i,0) + y_out(ii,2)=y_out2(i,2) + + y_out(ii,1)=y_out2(i,1)+y_out(i1,1) + y_out(ii,3)=y_out2(i,3)+y_out(i1,3) + y_out(ii,4)=y_out2(i,4)+y_out(i1,4) + + ii=ii+1 + ELSEIF(y_out2(i,0)>twopi)THEN + EXIT + ENDIF + ENDDO + len_y_out=ii-1 +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE patch_two_c +c----------------------------------------------------------------------- +c subprogram 23. yinterim. +c evaluates interpolated point at eta_i. assumes eta_i is within +c the eta-range of y_out +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE yinterim(y_out,len_y_out,eta_i,yinterim1) + + REAL(r8), DIMENSION(0:,0:), INTENT(IN) :: y_out + INTEGER, INTENT(IN) :: len_y_out + REAL(r8), INTENT(IN) :: eta_i + REAL(r8), DIMENSION(0:4), INTENT(OUT) :: yinterim1 + INTEGER :: i,j,eta_ii,k + TYPE(spline_type) :: yi + LOGICAL :: r_avg=.TRUE. +c----------------------------------------------------------------------- +c checking eta_i in range of y_out. +c----------------------------------------------------------------------- + IF(eta_iMAXVAL(y_out(0:len_y_out,0)))THEN + CALL program_stop("improper use of yinterim") + ENDIF +c----------------------------------------------------------------------- +c interpolating the numerical integral around an eta=0.0 point. +c I assume eta=0.0 won't be within the xpt_brackets since most +c tokamaks wont have an xpoint at the top or bottom of the LCFS +c----------------------------------------------------------------------- + IF(y_out(len_y_out,0)>eta_i)THEN + DO i=0,len_y_out,+1 + IF(y_out(i,0)>eta_i)THEN + eta_ii=i + EXIT + ENDIF + ENDDO + ELSEIF(y_out(len_y_out,0)eta1 +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE analytic_y_out(psifac,eta1,r1,eta2,r2,nstep2,yi,ydiffs, + $ outmat2f) + + REAL(r8), INTENT(IN) :: psifac,eta1,r1,eta2,r2 + INTEGER, INTENT(IN) :: nstep2 + REAL(r8), DIMENSION(0:,0:), INTENT(INOUT) :: yi + REAL(r8), DIMENSION(1:4,1:2,0:1) :: outmat + REAL(r8), DIMENSION(1:8) :: outmat2 + REAL(r8), DIMENSION(1:17), INTENT(OUT) :: outmat2f + LOGICAL :: debugL=.FALSE.,out=.FALSE. + REAL(r8), DIMENSION(1:2):: nu, xpt_dists + REAL(r8) :: r,z,eta,eval_BpOnBt,eval_dy1,eval_xpt_tol + REAL(r8), DIMENSION(0:nstepd,0:4) :: y_out1 + REAL(r8), DIMENSION(1:4), INTENT(OUT) :: ydiffs + INTEGER :: x_i,i,len_y1_out + TYPE(direct_bfield_type) :: bf + + yi=0.0 + ydiffs=0.0 + IF(eta2twopi/100 .OR. + $ ABS(xpt_gammas2(x_i)-xpt_gammas(x_i))>twopi/100) .AND. + $ (.NOT.outside_sep(x_i)))THEN + PRINT "(A)", "Straight x-point leg assumption bad for Xpt" + !PRINT "(i6)", x_i + !PRINT "(es16.10)", xpt_gammas(x_i)/pi + !PRINT "(es16.10)", xpt_gammas2(x_i)/pi + IF(.NOT.plot_convergence)THEN + CALL program_stop("fiddle with x-point settings.") + ENDIF + ENDIF +c----------------------------------------------------------------------- +c Setting initial values of yi. +c----------------------------------------------------------------------- + yi(0,0)=eta1 + yi(0,1)=0.0 + yi(0,2)=r1 + yi(0,3)=0.0 + yi(0,4)=0.0 +c----------------------------------------------------------------------- +c Print outputs directly if onecase is True: +c----------------------------------------------------------------------- +412 FORMAT(e20.12,",",e20.12,",",e20.12,",",e20.12,",",e20.12 + $",",e20.12,",",e20.12,",",e20.12,",",e20.12) + IF(onecase)CALL ascii_open(out_xpt_unit, + $ "xpt_tests/local_comp.csv","UNKNOWN") +c----------------------------------------------------------------------- +c Loop through eta values, computing analytic integrals +c----------------------------------------------------------------------- + DO i=1,nstep2,+1 + eta=eta1+(eta2-eta1)*(one*i/nstep2) + yi(i,0)=eta + + CALL direct_analytic_ints(x_i,r1,eta1 + $ ,eta,yi(i,1),yi(i,2),yi(i,3),yi(i,4),outmat,outmat2,debugL) + IF(i==nstep2)THEN + outmat2f(1:8)=outmat2(1:8) + ENDIF + + IF(.TRUE.)CALL direct_fl_int(psifac,xpt_brackets(1,1),eta, + $y_out1,len_y1_out,eval_BpOnBt,eval_dy1,eval_xpt_tol,.FALSE.,bf, + $ .TRUE.) + + WRITE(out_xpt_unit,412) + $ eta, + $ yi(i,1), + $ yi(i,2), + $ yi(i,3), + $ yi(i,4), + $ y_out1(len_y1_out,1), + $ y_out1(len_y1_out,2), + $ y_out1(len_y1_out,3), + $ y_out1(len_y1_out,4) + ENDDO + ydiffs(1)=(yi(nstep2,1)-y_out1(len_y1_out,1)) + $ /y_out1(len_y1_out,1) + ydiffs(2)=(yi(nstep2,2)-y_out1(len_y1_out,2)) + $ /y_out1(len_y1_out,2) + ydiffs(3)=(yi(nstep2,3)-y_out1(len_y1_out,3)) + $ /y_out1(len_y1_out,3) + ydiffs(4)=(yi(nstep2,4)-y_out1(len_y1_out,4)) + $ /y_out1(len_y1_out,4) + + outmat2f(9)=y_out1(len_y1_out,4) + outmat2f(10)=REAL(len_y1_out,8) !Number of integration steps.... + outmat2f(11)=y_out1(len_y1_out,1) + outmat2f(12)=y_out1(len_y1_out,2) + outmat2f(13)=y_out1(len_y1_out,3) + outmat2f(14)=yi(nstep2,1) + outmat2f(15)=yi(nstep2,2) + outmat2f(16)=yi(nstep2,3) + outmat2f(17)=yi(nstep2,4) + IF(onecase)CALL ascii_close(out_xpt_unit) +c----------------------------------------------------------------------- +c check how much r has shifted. +c----------------------------------------------------------------------- + IF(ABS(yi(nstep2,2)-r2)>r_tol)THEN + IF(.NOT. plot_convergence)WRITE(*,*) + $ "Analytic ints introduce radial discontinuity:" + IF(.NOT. plot_convergence)WRITE(*,*)ABS(yi(nstep2,2)-r2) + IF(.NOT. plot_convergence)CALL program_stop( + $ "Analytic integrals introduce a radial discontinuity > r_tol.") + ENDIF +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE analytic_y_out +c----------------------------------------------------------------------- +c subprogram 25. direct_Blocal. +c calculates local B-field displaced from some r,z point in polar +c coordinates. either Bnu or Brho depending on Bcase +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE direct_Blocal(r,z,nu,rho,Bcase,Bout) + + REAL(r8), INTENT(IN) :: r,z,nu,rho + CHARACTER, INTENT(IN) :: Bcase + REAL(r8), INTENT(OUT) :: Bout + REAL(r8) :: cosfac,sinfac + TYPE(direct_bfield_type) :: bf +c----------------------------------------------------------------------- +c calculation of Bout. +c----------------------------------------------------------------------- + cosfac=COS(nu) + sinfac=SIN(nu) + CALL direct_get_bfield(r+rho*cosfac,z+rho*sinfac,bf,1) + SELECT CASE (Bcase) + CASE ('n') + Bout=-sinfac*bf%br+cosfac*bf%bz + CASE ('r') + Bout=cosfac*bf%br+sinfac*bf%bz + CASE DEFAULT + CALL program_stop("Invalid Bcase in direct_Blocal") + END SELECT +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE direct_Blocal +c----------------------------------------------------------------------- +c subprogram 26. plot_xpt_convergence. +c calculates numerical vs analytic + numerical y_out integrals +c for a range of xpt_tol & BpBt_tol values, before printing outputs. +c used to generate the data that is plotted in +c https://doi.org/10.1088/1361-6587/add9ca Figs. 7 & 9. +c works on an equilibrium with 1 x-point. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE plot_xpt_convergence + + REAL(r8), DIMENSION(0:(4*nstepd+2*nstep2+6),0:4) :: y_outA,y_outI, + $ y_outDIFFS + REAL(r8), DIMENSION(0:(nstepd+2),0:4) :: y_out1 + REAL(r8) :: eval_BpOnBt,eval_dy1,eval_xpt_tol + TYPE(direct_bfield_type) :: bf + INTEGER :: len_y_out1,len_y_outA,len_y_outI + INTEGER :: i,j,ii + REAL(r8) :: psifac + + !Independent variables actual values + REAL(r8), DIMENSION(0:40,0:40,1:2) :: BpBt_evals,xpt_tol_evals, + $ dy1_evals + !Dependent variables + REAL(r8), DIMENSION(btsc1:btsc2,psfc1:psfc2,1:4) :: y_out_convgs + REAL(r8), DIMENSION(btsc1:btsc2,psfc1:psfc2,0:4) :: ff_convgs + REAL(r8), DIMENSION(btsc1:btsc2,psfc1:psfc2,1:17) :: outmat2big + REAL(r8), DIMENSION(1:4) :: difs + REAL(r8), DIMENSION(1:17) :: outmat2f + INTEGER, DIMENSION(2) :: x_starts + REAL(r8) :: psifacone,BpBt_tolone + REAL(r8), DIMENSION(btsc1:btsc2) :: BpBt_scan + REAL(r8), DIMENSION(psfc1:psfc2) :: psifac_vec + DO i=btsc1,btsc2,+1 !48 + BpBt_scan(i)=10.d0**(-i*0.1) + ENDDO + DO i=psfc1,psfc2,+1 !48 + !psifac_vec(i-8)=1.d0-10.d0**(-i*0.25) + psifac_vec(i)=1.d0-10.d0**(-i*0.2) + ENDDO + + psifacone=0.9987654321 + BpBt_tolone=0.012345 +c----------------------------------------------------------------------- +c initialisations: +c----------------------------------------------------------------------- + WRITE(*,*)"running convergence" + IF(num_xpts/=1)CALL program_stop("Need 1 x-point for convg. stud") + BpBt_evals=0.0 + xpt_tol_evals=0.0 + dy1_evals=0.0 + y_out_convgs=0.0 + ff_convgs=0.0 + + !xpt_tol=100 !will never be restained by angle + + DO i=btsc1,btsc2,+1 + WRITE(*,*)"BpBt_scan NEW:",BpBt_scan(i) + DO j=psfc1,psfc2,+1 + WRITE(*,*)i,j,"L0" + psifac=psifac_vec(j) + WRITE(*,*)"Psifac, BpBt:" + WRITE(*,*)psifac,BpBt_scan(i) +c----------------------------------------------------------------------- +c GENERATING THE NUMERICALLY INTEGRATED y_out:: +c----------------------------------------------------------------------- + IF(onecase)THEN + psifac=psifacone + ENDIF + etol=1.d-14 !lowest we can go? + nstepd=20000 !big KEEP UNDER 20000 + BpBt_tol=1d-99 !we will never 'snag' an x-point' + CALL direct_fl_int(psifac,zero,twopi,y_outI,len_y_outI + $ ,eval_BpOnBt,eval_dy1,eval_xpt_tol,.FALSE.,bf,.FALSE.) +c----------------------------------------------------------------------- +c GENERATING THE ANALYTICALLY INTEGRATED y_out:: +c ASSUME 1xpt equilibrium +c----------------------------------------------------------------------- + BpBt_tol=BpBt_scan(i) + IF(onecase)THEN + BpBt_tol=BpBt_tolone + ENDIF + CALL direct_fl_int(psifac,zero,twopi,y_out1,len_y_out1 + $ ,eval_BpOnBt,eval_dy1,eval_xpt_tol,.TRUE.,bf,.FALSE.) +c----------------------------------------------------------------------- +c GENERATING THE COMBINED y_out (analytic + numerical):: +c----------------------------------------------------------------------- + IF(y_out1(len_y_out1,0)