module coefficients use constants, only: dp, pi, machine_eps use fieldline_mod, only: flock_of_fieldlines_t use, intrinsic :: ieee_arithmetic, only: ieee_is_nan implicit none private public :: calc_finite_boundary_layer_correction public :: calc_gradient_scaling_factor_r_eff public :: calc_offset_coefficients, calc_nu_star_crit contains !> Finite boundary layer correction Lambda_S to the bootstrap current. !! !! Accounts for the finite width of the trapped-passing boundary layer !! (Helander, Geiger and Maaßberg, Phys. Plasmas 2011) function calc_finite_boundary_layer_correction(flock, & field, & R, & dr_dAtheta) & result(Lambda_S) use surface_average_mod, only: surface_average_t, calc_surface_averages use field_base, only: field_t type(flock_of_fieldlines_t), intent(in) :: flock class(field_t), intent(in) :: field real(dp), intent(in) :: R !! major radius in metres real(dp), intent(in) :: dr_dAtheta !! e.g. from calc_gradient_scaling_factor_r_eff real(dp) :: Lambda_S real(dp) :: I_ref, M_pol, N_tor, iota, eta_b real(dp) :: B_theta_covariant, B_phi_covariant real(dp) :: I_ref_hat real(dp) :: helical_factor type(surface_average_t) :: average if (R < machine_eps) then print *, "error: major radius R must be positive." print *, "R =", R error stop end if if (ieee_is_nan(dr_dAtheta)) then print *, "error: dr_dAtheta is NaN." error stop end if call field%get_covariant_components(B_theta_covariant, B_phi_covariant) if (ieee_is_nan(B_theta_covariant)) then print *, "error: B_theta_covariant is NaN." error stop end if if (ieee_is_nan(B_phi_covariant)) then print *, "error: B_phi_covariant is NaN." error stop end if call calc_surface_averages(flock, average) M_pol = flock%M_pol N_tor = flock%N_tor iota = flock%iota I_ref = flock%I_ref eta_b = flock%eta_b helical_factor = (B_phi_covariant*M_pol + B_theta_covariant*N_tor) & /(M_pol*iota - N_tor) I_ref_hat = I_ref/eta_b/(2.0_dp*pi*R) Lambda_S = sqrt(I_ref_hat)*average%B_squared*eta_b**2/average%lambda_b Lambda_S = 1.5_dp*0.855_dp/sqrt(pi)*Lambda_S Lambda_S = -helical_factor*dr_dAtheta*Lambda_S if (ieee_is_nan(Lambda_S)) then print *, "error: Lambda_S is NaN." error stop end if end function calc_finite_boundary_layer_correction !> Scaling factor dr/dA_theta for choice of effective radius !! !! Computes the factor to convert bootstrap coefficients so that !! they can be used with gradients of the "effective radius" r, defined by !! dV/dr = S, where V is the enclosed volume and S the flux surface area. function calc_gradient_scaling_factor_r_eff(flock, psi_tor_edge, sign_sqrtg) & result(dr_dAtheta) use surface_average_mod, only: surface_average_t, calc_surface_averages type(flock_of_fieldlines_t), intent(in) :: flock real(dp), intent(in) :: psi_tor_edge !! total toroidal flux at the plasma edge in Weber integer, intent(in) :: sign_sqrtg !! sign of the Boozer Jacobian: +1 or -1 real(dp) :: dr_dAtheta type(surface_average_t) :: average if (abs(sign_sqrtg) /= 1) then print *, "error: sign_sqrtg must be +1 or -1." print *, "sign_sqrtg =", sign_sqrtg error stop end if call calc_surface_averages(flock, average) if (abs(average%nabla_s*psi_tor_edge) < machine_eps) then print *, "error: nabla_s * psi_tor_edge is zero." print *, "nabla_s =", average%nabla_s print *, "psi_tor_edge =", psi_tor_edge error stop end if dr_dAtheta = real(sign_sqrtg, dp)/(average%nabla_s*psi_tor_edge) if (ieee_is_nan(dr_dAtheta)) then print *, "error: dr_dAtheta is NaN." error stop end if end function calc_gradient_scaling_factor_r_eff !> Compute offset coefficients Lambda_A and Lambda_B. !! !! Coefficients are internally computed for gradients in respect to the !! poloidal Boozer vector-potential component i.e. !! boostrap current = coefficient * dp/dA_theta !! One needs to apply a conversion factor dr_dA_theta to use them with !! gradients of other surface labels r. Furthermore, they are given for a !! specific collisionality defintion which is set by major radius R. subroutine calc_offset_coefficients(flock, R, dr_dAtheta, & Lambda_A, Lambda_B) use deviation, only: calc_deviation use field_base, only: field_t type(flock_of_fieldlines_t), intent(in) :: flock real(dp), intent(in) :: R !! major radius in metres real(dp), intent(in) :: dr_dAtheta !! defining the surface label r real(dp), intent(out) :: Lambda_A, Lambda_B real(dp) :: deviation_A, deviation_B call calc_deviation(flock, deviation_A, deviation_B) Lambda_A = deviation_A*dr_dAtheta*sqrt(0.5_dp*R*pi) Lambda_B = deviation_B*0.5_dp*R*pi*dr_dAtheta end subroutine calc_offset_coefficients !> Critical collisionality marking the lower limit of validity for the !! asymptotic model. Collisionality is defined by major radius R. function calc_nu_star_crit(flock, R) result(nu_star_crit) type(flock_of_fieldlines_t), intent(in) :: flock real(dp), intent(in) :: R !! major radius in metres real(dp) :: nu_star_crit real(dp) :: eta_b, I_ref_hat, max_delta_eta if (R < machine_eps) then print *, "error: major radius R must be positive." print *, "R =", R error stop end if eta_b = flock%eta_b I_ref_hat = flock%I_ref/eta_b/(2.0_dp*pi*R) max_delta_eta = eta_b - 1.0_dp/minval(flock%fieldlines%B_max(1)) nu_star_crit = 0.125_dp*(max_delta_eta/eta_b)**2.0_dp/I_ref_hat if (ieee_is_nan(nu_star_crit)) then print *, "error: Lambda_S is NaN." error stop end if end function calc_nu_star_crit end module coefficients