module safety_factor_m !! Computes safety factor !! N.B. This module is currently axisymmetric - can only calculate !! pitch and safety factor for 2D equilibria !! PHITEMP parameter needed for equilibrium interface use precision_m, only : FP, DP, DP_EPS use constants_m, only : TWO_PI use screen_io_m, only: get_stderr use comm_handling_m, only: is_master use equilibrium_m, only: equilibrium_t use coords_polar_m, only : polar_to_cart use dop853_module use dop853_constants implicit none private :: set_equilibrium private :: pitch public :: safety_factor type, private, extends(dop853_class) :: dop853_pitch_t !! Tracer for q-profile based on pitch angle class(equilibrium_t), private, pointer :: equi => null() real(FP), public :: rho contains procedure, public :: set_equilibrium procedure, public :: fcn => pitch procedure, public :: solout => solout_none end type contains subroutine set_equilibrium(self, equi) !! Sets equilibrium to ode integrator class(dop853_pitch_t), intent(inout) :: self class(equilibrium_t), target, intent(in) :: equi self%equi => equi end subroutine subroutine pitch(me, x, y, f) !! Computes B^phi/B^theta at position theta, used as integrand for dop853 class(dop853_pitch_t),intent(inout) :: me real(DP), intent(in) :: x !! Poloidal angle theta real(DP), dimension(:), intent(in) :: y !! Current value of integral real(DP), dimension(:), intent(out) :: f !! Integrand = pitch angle real(FP) :: theta, xc, yc, rtorus, jac, pitch_fp real(FP), parameter :: PHITEMP = 0.0_FP !! pitch interfaces with dop853, thus phi can't be regular intent(in). !! This routine needs further updating to be 3D theta = real(x, kind=FP) call polar_to_cart(me%equi, me%rho, theta, PHITEMP, xc, yc) rtorus = sqrt( (xc - me%equi%x0)**2 + (yc - me%equi%y0)**2 ) jac = me%equi%jacobian(xc, yc, PHITEMP) pitch_fp = me%equi%btor(xc, yc, PHITEMP) / sqrt(me%equi%bx(xc, yc, PHITEMP)**2 + & me%equi%by(xc, yc, PHITEMP)**2)* rtorus / jac f(1) = real(pitch_fp, kind=DP) end subroutine subroutine solout_none(me, nr, xold, x, y, irtrn, xout) !! An empty function needed in order to run the field line tracer class(dop853_pitch_t),intent(inout) :: me !! Instance of the type integer,intent(in) :: nr !! Grid point (0,1,...) real(DP),intent(in) :: xold !! The preceding grid point real(DP),intent(in) :: x !! Current grid point real(DP),dimension(:),intent(in) :: y !! State vector \( y(x) \) [size n] integer,intent(inout) :: irtrn !! Serves to interrupt the integration, if irtrn is set <0 real(DP),intent(out) :: xout !! Can be used for efficient intermediate output end subroutine real(FP) function safety_factor(equi, rho) !! Compute safety factor q(rho) class(equilibrium_t), target, intent(inout) :: equi !! Equilibrium defining the field line real(FP), intent(in) :: rho !! Flux surface label type(dop853_pitch_t) :: odeint logical :: odestat_ok integer :: idid real(DP), dimension(1) :: rtol, atol real(DP) :: t, tf real(DP), dimension(1) :: qsafe call odeint%set_equilibrium(equi) odeint%rho = rho ! q profile restricted by absoulte error) rtol(1) = DP_EPS atol(1) = 1.0E-8_DP call odeint%initialize(n = 1, & iprint = get_stderr(), & status_ok = odestat_ok ) if (.not.odestat_ok) then write(get_stderr(), *) & 'error(saftey_factor) :: ode initialisation' error stop endif t = 0.0_FP tf = TWO_PI qsafe(1) = 0.0_FP call odeint%integrate(t, qsafe, tf, rtol, atol, iout=0, idid=idid) if (idid /= 1) then write(get_stderr(), *) & 'error(saftey_factor) :: dop853 did not converge, idid=', idid error stop endif safety_factor = real(qsafe(1), kind=FP) / TWO_PI end function end module