module gyroviscosity_m !! Computations related with gyroviscosity use precision_grillix_m, only : GP, GP_EPS use error_handling_grillix_m, only: handle_error use status_codes_grillix_m, only : GRILLIX_ERR_OTHER use descriptors_m, only : BND_TYPE_NEUMANN use comm_handler_m, only : comm_handler_t use equilibrium_m, only : equilibrium_t use equilibrium_storage_m, only : equilibrium_storage_t use boundaries_perp_m, only : set_perpbnds, extrapolate_ghost_points use mesh_cart_m, only : mesh_cart_t use parallel_map_m, only : parallel_map_t use variable_m, only: variable_t use inplane_operators_m, only : inplane_operators_t use params_brag_model_m, only : eta_i0, tratio, delta, beta, heat_fac use params_brag_pardiss_model_m, only : free_streaming_limiter_qfac, aspect_ratio_inv use params_brag_switches_m, only : swb_neocl_visc, swb_flutter_visc, swb_upar_flutter_viscion_grad implicit none type, public :: gyroviscosity_t !! Holds information and performs computations related with gyroviscosity type(variable_t), private :: uparb32 !! Parallel velocity time B^(3/2) type(variable_t), private :: qiparb32 !! Parallel conductive ion heat flux times B^(3/2) real(GP), allocatable, dimension(:), public :: eta_flow_v !! Values of flow viscosity (on canonical mesh) real(GP), allocatable, dimension(:), private :: eta_heat_v !! Values of heat viscosity (on canonical mesh) contains procedure, public :: init => init_gyroviscosity procedure, public :: drive => drive_gyroviscosity procedure, public :: eval_gstress procedure, public :: eval_curvature procedure, public :: eval_pargrad_b32 final :: destructor end type public :: eta_neo_flow private :: eta_neo_heat private :: neo_coll_time contains subroutine init_gyroviscosity(self, comm_handler, mesh_cano, mesh_stag) !! Initialise gyroviscosity class(gyroviscosity_t), intent(inout) :: self !! Instance of type type(comm_handler_t), intent(in) :: comm_handler !! Communicators type(mesh_cart_t), intent(in) :: mesh_cano !! Mesh (canonical) within poloidal plane type(mesh_cart_t), intent(in) :: mesh_stag !! Mesh (staggered) within poloidal plane ! Initialse variable, needed to be communicated call self%uparb32%init(comm_handler, mesh_stag%get_n_points(), staggered=.true.) call self%uparb32%create_halos(comm_handler, 0, 1) call self%qiparb32%init(comm_handler, mesh_stag%get_n_points(), staggered=.true.) call self%qiparb32%create_halos(comm_handler, 0, 1) allocate(self%eta_flow_v(mesh_cano%get_n_points())) allocate(self%eta_heat_v(mesh_cano%get_n_points())) ! Check consistency of input parameters if ((swb_neocl_visc > 0.0_GP) .and. (aspect_ratio_inv < GP_EPS)) then call handle_error('Running with neoclassical corrections (swb_neocl_visc>0), & aspect_ratio_inv must be set larger than zero', & GRILLIX_ERR_OTHER, __LINE__, __FILE__) endif end subroutine subroutine drive_gyroviscosity(self, comm_handler, mesh_cano, mesh_stag, equi_on_stag, & ne_v, ti_v, upar_v, qipar_v) !! Performs auxiliary precomputations necessary to evaluate guroviscosity !! - Precomputes viscosities !! - Communicates u*B^3/2 and qipar*B^(3/2) class(gyroviscosity_t), intent(inout) :: self !! Instance of type type(comm_handler_t), intent(in) :: comm_handler !! Communicators type(mesh_cart_t), intent(in) :: mesh_cano !! Mesh (canonical) within poloidal plane type(mesh_cart_t), intent(in) :: mesh_stag !! Mesh (staggered) within poloidal plane class(equilibrium_storage_t), intent(in) :: equi_on_stag !! Equilibrim quantities on staggered mesh real(GP), dimension(mesh_cano%get_n_points()), intent(in) :: ne_v !! Values of density real(GP), dimension(mesh_cano%get_n_points()), intent(in) :: ti_v !! Values of ion temperature real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: upar_v !! Values of parallel ion velocity real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: qipar_v !! Values of parallel ion heat flux, for comoputation of heat-flux related viscosity integer :: l !$omp parallel default(none) private(l) & !$omp shared(self, mesh_cano, mesh_stag, equi_on_stag, ne_v, ti_v, upar_v, qipar_v) !$omp do do l = 1, mesh_cano%get_n_points() self%eta_flow_v(l) = eta_neo_flow(ne_v(l), ti_v(l)) self%eta_heat_v(l) = eta_neo_heat(ne_v(l), ti_v(l)) enddo !$omp end do !$omp do do l = 1, mesh_stag%get_n_points() self%uparb32%vals(l) = upar_v(l) * sqrt(equi_on_stag%absb(l)**3) self%qiparb32%vals(l) = qipar_v(l) * sqrt(equi_on_stag%absb(l)**3) enddo !$omp end do !$omp end parallel call self%uparb32%start_comm_halos(comm_handler) call self%uparb32%finalize_comm_halos(comm_handler) call self%qiparb32%start_comm_halos(comm_handler) call self%qiparb32%finalize_comm_halos(comm_handler) end subroutine subroutine destructor(self) !! Destructor type(gyroviscosity_t), intent(inout) :: self !! Instance of type if (allocated(self%eta_flow_v)) deallocate(self%eta_flow_v) if (allocated(self%eta_heat_v)) deallocate(self%eta_heat_v) end subroutine subroutine eval_gstress(self, mesh_cano, equi_on_cano, opsinplane_cano, map, & ne_v, pot_v, ti_v, apar_fluct_cano_v, gstress, & eval_flow_crv, eval_flow_par, eval_heat_crv, eval_heat_par) !! Computation of ion stress function G = G_flow + G_heat !! with "beyond Braginskii extensions" according to parameters class(gyroviscosity_t), intent(in) :: self !! Instance of type type(mesh_cart_t), intent(in) :: mesh_cano !! Mesh (canonical) within poloidal plane class(equilibrium_storage_t), intent(in) :: equi_on_cano !! Equilibrim quantities on canonical mesh type(inplane_operators_t), intent(in) :: opsinplane_cano !! In-plane operators (canonical) type(parallel_map_t), intent(in) :: map !! Parallel map real(GP), dimension(mesh_cano%get_n_points()), intent(in) :: ne_v !! Electron density real(GP), dimension(mesh_cano%get_n_points()), intent(in) :: pot_v !! Electrostatic potential real(GP), dimension(mesh_cano%get_n_points()), intent(in) :: ti_v !! Ion temperature real(GP), dimension(mesh_cano%get_n_points()), intent(in) :: apar_fluct_cano_v !! Fluctuation of apar used for flutter operators !! It is assumed that communication has already been done to do mapping to canonical mesh real(GP), dimension(mesh_cano%get_n_points()), intent(out) :: gstress !! Values of ion viscous stress function G logical, intent(in), optional :: eval_flow_crv !! If true computes part of G_flow arising from curvature terms !! (default = T) !! This feature is useful for technical reasons, !! when only evaluation of individual parts of G is needed logical, intent(in), optional :: eval_flow_par !! If true computes computes part of G_flow arising from parallel divergence of upar !! (default = T) logical, intent(in), optional :: eval_heat_crv !! If true computes part of G_heat arising from curvature terms !! (default = T) logical, intent(in), optional :: eval_heat_par !! If true computes part of G_flow arising from paralle divergence of qipar !! (default = T) logical :: apply_flow_crv, apply_flow_par, apply_heat_crv, apply_heat_par integer :: l, ki, kb real(GP) :: crv_pot, crv_ti, crv_pi, & duparb32dpar, duparb32dpar_flutter, dqiparb32dpar, dqiparb32dpar_flutter integer, dimension(mesh_cano%get_n_points_boundary()) :: bnd_descrs_nmn real(GP), dimension(mesh_cano%get_n_points()) :: uparb32_cano, qiparb32_cano ! Select, which parts to compute apply_flow_crv = .true. if (present(eval_flow_crv)) apply_flow_crv = eval_flow_crv apply_flow_par = .true. if (present(eval_flow_par)) apply_flow_par = eval_flow_par apply_heat_crv = .true. if (present(eval_heat_crv)) apply_heat_crv = eval_heat_crv apply_heat_par = .true. if (present(eval_heat_par)) apply_heat_par = eval_heat_par if ((heat_fac < GP_EPS) .or. & (swb_neocl_visc < GP_EPS)) then apply_heat_crv = .false. apply_heat_par = .false. endif !$omp parallel default(none) & !$omp private(ki, kb, l, crv_pot, crv_pi, crv_ti, & !$omp duparb32dpar, duparb32dpar_flutter, dqiparb32dpar, dqiparb32dpar_flutter) & !$omp shared(self, mesh_cano, equi_on_cano, map, opsinplane_cano, swb_flutter_visc, & !$omp tratio, delta, beta, bnd_descrs_nmn, & !$omp apply_flow_crv, apply_flow_par, apply_heat_crv, apply_heat_par, & !$omp ne_v, pot_v, ti_v, uparb32_cano, qiparb32_cano, & !$omp apar_fluct_cano_v, gstress) if (apply_flow_par) then call map%flat_stag_to_cano(self%uparb32%vals, self%uparb32%hbwd, uparb32_cano) endif if(apply_heat_par) then call map%flat_stag_to_cano(self%qiparb32%vals, self%qiparb32%hbwd, qiparb32_cano) endif !$omp do do ki = 1, mesh_cano%get_n_points_inner() l = mesh_cano%inner_indices(ki) gstress(l) = 0.0_GP if (apply_flow_crv) then crv_pot = opsinplane_cano%curvature(pot_v, l) crv_ti = opsinplane_cano%curvature(ti_v, l) crv_pi = ti_v(l) * opsinplane_cano%curvature(ne_v, l) & + ne_v(l) * crv_ti gstress(l) = gstress(l) + 0.5_GP * self%eta_flow_v(l) * (crv_pot + tratio * crv_pi / ne_v(l)) endif if (apply_heat_crv) then gstress(l) = gstress(l) + 1.25_GP * self%eta_heat_v(l) * tratio * crv_ti endif if (apply_flow_par) then duparb32dpar = map%par_divb_single(self%uparb32%vals, self%uparb32%hbwd, l) duparb32dpar_flutter = swb_flutter_visc & * opsinplane_cano%flutter_div(apar_fluct_cano_v, uparb32_cano, l) gstress(l) = gstress(l) - self%eta_flow_v(l) * 2.0_GP * (duparb32dpar + duparb32dpar_flutter) & / sqrt(equi_on_cano%absb(l)**3) endif if (apply_heat_par) then dqiparb32dpar = map%par_divb_single(self%qiparb32%vals, self%qiparb32%hbwd, l) dqiparb32dpar_flutter = swb_flutter_visc & * opsinplane_cano%flutter_div(apar_fluct_cano_v, qiparb32_cano, l) gstress(l) = gstress(l) - self%eta_heat_v(l) * 2.0_GP * (dqiparb32dpar + dqiparb32dpar_flutter) & / (sqrt(equi_on_cano%absb(l)**3) * ne_v(l) * ti_v(l)) endif enddo !$omp end do ! Extrapolate G to boundary and ghost points ! Note: Could this be handled better, especially w.r.t MMS !$omp do do kb = 1, mesh_cano%get_n_points_boundary() bnd_descrs_nmn(kb) = BND_TYPE_NEUMANN enddo !$omp end do call set_perpbnds(mesh_cano, bnd_descrs_nmn, gstress) call extrapolate_ghost_points(mesh_cano, gstress) !$omp end parallel end subroutine subroutine eval_curvature(self, mesh_cano, equi_on_cano, opsinplane_cano, map, & ne_v, pot_v, ti_v, apar_fluct_cano_v, crv_gstress) !! Evaluates curvature of ion stress function. !! This routine is present, because there is an issue if one simply evaluated the curvature !! of gstress obtained via ion_stress_function. !! Squared curvature operators appear that would imply chequerboarding. !! Therefore the evaluation of squared curvature operators has to be performed numerically !! careful and explicit. class(gyroviscosity_t), intent(in) :: self !! Instance of type type(mesh_cart_t), intent(in) :: mesh_cano !! Mesh (canonical) within poloidal plane class(equilibrium_storage_t), intent(in) :: equi_on_cano !! Equilibrim quantities on canonical mesh type(inplane_operators_t), intent(in) :: opsinplane_cano !! In-plane operators (canonical) type(parallel_map_t), intent(in) :: map !! Parallel map real(GP), dimension(mesh_cano%get_n_points()), intent(in) :: ne_v !! Electron density real(GP), dimension(mesh_cano%get_n_points()), intent(in) :: pot_v !! Electrostatic potential real(GP), dimension(mesh_cano%get_n_points()), intent(in) :: ti_v !! Ion temperature real(GP), dimension(mesh_cano%get_n_points()), intent(in) :: apar_fluct_cano_v !! Fluctuation of apar on canonical mesh used for flutter operators real(GP), dimension(mesh_cano%get_n_points_inner()), intent(out) :: crv_gstress !! Values of curvature of ion viscous stress function G logical :: eval_heat integer :: ki, l real(GP) :: crv_gpar, crv_eta_flow, crv_eta_heat, crv_pot, crv_ti, crv_ne, crv_pi, & crvsqrd_pot, crvsqrd_ne, crvsqrd_ti, crvsqrd_pi real(GP), dimension(mesh_cano%get_n_points()) :: gstress_par ! Determine if heat conduction part is computed if (swb_neocl_visc > 0.0_GP .and. & heat_fac > 0.0_GP) then eval_heat = .true. else eval_heat = .false. endif ! Compute parallel part of viscous stress tensor call self%eval_gstress(mesh_cano, equi_on_cano, opsinplane_cano, map, & ne_v, pot_v, ti_v, apar_fluct_cano_v, gstress_par, & eval_flow_crv=.false., eval_flow_par=.true., & eval_heat_crv=.false., eval_heat_par=eval_heat) !$omp parallel default(none) & !$omp private(ki, l, crv_gpar, crv_eta_flow, crv_eta_heat, crv_pot, crv_ti, & !$omp crv_ne, crv_pi, crvsqrd_pot, crvsqrd_ne, crvsqrd_ti, crvsqrd_pi) & !$omp shared(self, mesh_cano, opsinplane_cano, tratio, eval_heat, ne_v, ti_v, pot_v, & !$omp gstress_par, crv_gstress) !$omp do do ki = 1, mesh_cano%get_n_points_inner() l = mesh_cano%inner_indices(ki) ! Curvatures crv_gpar = opsinplane_cano%curvature(gstress_par, l) crv_eta_flow = opsinplane_cano%curvature(self%eta_flow_v, l) crv_pot = opsinplane_cano%curvature(pot_v, l) crv_ti = opsinplane_cano%curvature(ti_v, l) crv_ne = opsinplane_cano%curvature(ne_v, l) crv_pi = ti_v(l) * crv_ne + ne_v(l) * crv_ti crvsqrd_pot = opsinplane_cano%curvature_sqrd(pot_v, l) crvsqrd_ne = opsinplane_cano%curvature_sqrd(ne_v, l) crvsqrd_ti = opsinplane_cano%curvature_sqrd(ti_v, l) crvsqrd_pi = ne_v(l) * crvsqrd_ti + ti_v(l) * crvsqrd_ne & + 2.0_GP * crv_ne * crv_ti crv_gstress(ki) = crv_gpar & + 0.5_GP * crv_eta_flow * (crv_pot + tratio * crv_pi / ne_v(l)) & + 0.5_GP * self%eta_flow_v(l) * (crvsqrd_pot & + tratio * crvsqrd_pi / ne_v(l) & - tratio * crv_ne * crv_pi / ne_v(l)**2) if (eval_heat) then crv_eta_heat = opsinplane_cano%curvature(self%eta_heat_v, l) crv_gstress(ki) = crv_gstress(ki) & + tratio * 1.25_GP * (crv_eta_heat * crv_ti) & + tratio * 1.25_GP * self%eta_heat_v(l) * crvsqrd_ti endif end do !$omp end do !$omp end parallel end subroutine subroutine eval_pargrad_b32(self, comm_handler, mesh_cano, mesh_stag, equi_on_cano, equi_on_stag, & opsinplane_cano, opsinplane_stag, map, & ne_v, pot_v, ti_v, apar_fluct_v, apar_fluct_cano_v, pargrad_gstress, & eval_flow_crv, eval_flow_par, eval_heat_crv, eval_heat_par) !! Evaluates B^(3/2)par_grad(G/B^(3/2)) of ion stress function. !! Used for the parallel momentum equation class(gyroviscosity_t), intent(in) :: self !! Instance of type type(comm_handler_t), intent(in) :: comm_handler !! Communicators type(mesh_cart_t), intent(in) :: mesh_cano !! Mesh (canonical) within poloidal plane type(mesh_cart_t), intent(in) :: mesh_stag !! Mesh (staggered) within poloidal plane class(equilibrium_storage_t), intent(in) :: equi_on_cano !! Equilibrim quantities on canonical mesh class(equilibrium_storage_t), intent(in) :: equi_on_stag !! Equilibrim quantities on staggered mesh type(inplane_operators_t), intent(in) :: opsinplane_cano !! In-plane operators (canonical) type(inplane_operators_t), intent(in) :: opsinplane_stag !! In-plane operators (staggered) type(parallel_map_t), intent(in) :: map !! Parallel map real(GP), dimension(mesh_cano%get_n_points()), intent(in) :: ne_v !! Electron density real(GP), dimension(mesh_cano%get_n_points()), intent(in) :: pot_v !! Electrostatic potential real(GP), dimension(mesh_cano%get_n_points()), intent(in) :: ti_v !! Ion temperature real(GP), dimension(mesh_stag%get_n_points()), intent(in) :: apar_fluct_v !! Fluctuation of apar used for flutter operators real(GP), dimension(mesh_cano%get_n_points()), intent(in) :: apar_fluct_cano_v !! Fluctuation of apar on canonical mesh used for flutter operators real(GP), dimension(mesh_stag%get_n_points_inner()), intent(out) :: pargrad_gstress !! Values of ion viscous stress function G logical, intent(in) :: eval_flow_crv !! If true computes part of G_flow arising from curvature terms !! This feature is useful for technical reasons, !! when only evaluation of individual parts of G is needed logical, intent(in) :: eval_flow_par !! If true computes computes part of G_flow arising from parallel divergence of upar logical, intent(in) :: eval_heat_crv !! If true computes part of G_heat arising from curvature terms logical, intent(in) :: eval_heat_par !! If true computes computes part of G_flow arising from parallel divergence of qipar integer :: ki, l real(GP) :: dgstress32dpar, dgstress32dpar_flutter type(variable_t) :: gstress32 real(GP), dimension(mesh_stag%get_n_points()) :: gstress32_stag call gstress32%init(comm_handler, mesh_cano%get_n_points(), staggered=.false.) call gstress32%create_halos(comm_handler, 1, 0) call self%eval_gstress(mesh_cano, equi_on_cano, opsinplane_cano, map, & ne_v, pot_v, ti_v, apar_fluct_cano_v, gstress32%vals, & eval_flow_crv, eval_flow_par, eval_heat_crv, eval_heat_par) !$omp parallel default(none) private(l) shared(mesh_cano, equi_on_cano, gstress32) !$omp do do l = 1, mesh_cano%get_n_points() gstress32%vals(l) = gstress32%vals(l) / sqrt(equi_on_cano%absb(l)**3) enddo !$omp end do !$omp end parallel call gstress32%start_comm_halos(comm_handler) call gstress32%finalize_comm_halos(comm_handler) !$omp parallel default(none) & !$omp private(ki, l, dgstress32dpar, dgstress32dpar_flutter) & !$omp shared(mesh_stag, equi_on_stag, map, opsinplane_stag, & !$omp swb_upar_flutter_viscion_grad, delta, beta, & !$omp gstress32, gstress32_stag, apar_fluct_v, pargrad_gstress) call map%lift_cano_to_stag(gstress32%vals, gstress32%hfwd, gstress32_stag) !$omp do do ki = 1, mesh_stag%get_n_points_inner() l = mesh_stag%inner_indices(ki) dgstress32dpar = map%par_grad_single(gstress32%vals, gstress32%hfwd, l) dgstress32dpar_flutter = swb_upar_flutter_viscion_grad & * opsinplane_stag%flutter_grad(apar_fluct_v, gstress32_stag, l) pargrad_gstress(ki) = (dgstress32dpar + dgstress32dpar_flutter) * sqrt(equi_on_stag%absb(l)**3) enddo !$omp end do !$omp end parallel end subroutine pure function neo_coll_time(ne_l, ti_l) result(res) !! Computes neoclassical collison time real(GP), intent(in) :: ne_l !! Density value real(GP), intent(in) :: ti_l !! Ion temperature value real(GP) :: res real(GP) :: eps32 eps32 = sqrt(aspect_ratio_inv**3) res = sqrt(2.0_GP * tratio) * eps32 * eta_i0 * ti_l**2 & / (0.96_GP * free_streaming_limiter_qfac * ne_l) end function pure function eta_neo_flow(ne_l, ti_l) result(res) !! Computes ion flow viscosity !! with "beyond Braginskii extensions" according to parameters real(GP), intent(in) :: ne_l !! Density value real(GP), intent(in) :: ti_l !! Ion temperature value real(GP) :: res real(GP) :: eps32, tau_neo, fc res = eta_i0 * sqrt(ti_l**5) if (swb_neocl_visc > GP_EPS) then eps32 = sqrt(aspect_ratio_inv**3) tau_neo = neo_coll_time(ne_l, ti_l) fc = 1.0_GP / (1.0_GP + swb_neocl_visc * tau_neo / eps32) & * 1.0_GP / (1.0_GP + swb_neocl_visc * tau_neo) res = res * fc endif end function pure function eta_neo_heat(ne_l, ti_l) result(res) !! Computes ion heat viscosity !! "beyond Braginskii extensions" real(GP), intent(in) :: ne_l !! Density value real(GP), intent(in) :: ti_l !! Ion temperature value real(GP) :: res real(GP), parameter :: rt2 = sqrt(2.0_GP) real(GP) :: tau_neo, nurt, nusq, eps3, eta_flow, kt if ((swb_neocl_visc < GP_EPS) .or. (heat_fac < GP_EPS)) then res = 0.0_GP return endif tau_neo = neo_coll_time(ne_l, ti_l) nurt = sqrt(rt2 / tau_neo) nusq = (rt2 / tau_neo)**2 eps3 = aspect_ratio_inv**3 eta_flow = eta_neo_flow(ne_l, ti_l) kt = ( -0.17_GP + 1.05_GP * nurt + & 2.7_GP * nusq * eps3) / & ( 1.0_GP + 0.7_GP * nurt + nusq * eps3 ) res = swb_neocl_visc * heat_fac * 8.0_GP / 15_GP * (kt - 1.0_GP) * eta_flow end function end module