module boundaries_braginskii_m ! Stores information on boundary conditions of fields for Braginskii model use polars_m, only : polars_t use equilibrium_m, only : equilibrium_t use comm_handler_m, only : comm_handler_t use precision_grillix_m, only : GP, GP_NAN use error_handling_grillix_m, only: handle_error, error_info_t use status_codes_grillix_m, only : GRILLIX_ERR_OTHER use descriptors_braginskii_m, only : BND_BRAGTYPE_ZONAL_NEUMANN ! From PARALLAX use screen_io_m, only : get_stdout use mesh_cart_m, only : mesh_cart_t use descriptors_m, only : BND_TYPE_DIRICHLET_ZERO, BND_TYPE_NEUMANN, & DISTRICT_CORE, DISTRICT_WALL, DISTRICT_DOME, DISTRICT_OUT use interpolation_m, only: interpol1d use polar_operators_m, only : surface_average use boundaries_perp_m, only : set_perpbnds !! Parameters use params_brag_boundaries_perp_m, only : & bnddescr_ne_core, bnddescr_ne_wall, bnddescr_ne_dome, bnddescr_ne_out, & bnddescr_pot_core, bnddescr_pot_wall, bnddescr_pot_dome, bnddescr_pot_out, & bnddescr_vort_core, bnddescr_vort_wall, bnddescr_vort_dome, bnddescr_vort_out, & bnddescr_apar_core, bnddescr_apar_wall, bnddescr_apar_dome, bnddescr_apar_out, & bnddescr_jpar_core, bnddescr_jpar_wall, bnddescr_jpar_dome, bnddescr_jpar_out, & bnddescr_te_core, bnddescr_te_wall, bnddescr_te_dome, bnddescr_te_out, & bnddescr_ti_core, bnddescr_ti_wall, bnddescr_ti_dome, bnddescr_ti_out, & bnddescr_upar_core, bnddescr_upar_wall, bnddescr_upar_dome, bnddescr_upar_out, & bnddescr_qe_core, bnddescr_qe_wall, bnddescr_qe_dome, bnddescr_qe_out, & bnddescr_qi_core, bnddescr_qi_wall, bnddescr_qi_dome, bnddescr_qi_out, & bndval_ne_core, bndval_ne_wall, bndval_ne_dome, bndval_ne_out, & bndval_pot_core, bndval_pot_wall, bndval_pot_dome, bndval_pot_out, & bndval_vort_core, bndval_vort_wall, bndval_vort_dome, bndval_vort_out, & bndval_apar_core, bndval_apar_wall, bndval_apar_dome, bndval_apar_out, & bndval_jpar_core, bndval_jpar_wall, bndval_jpar_dome, bndval_jpar_out, & bndval_te_core, bndval_te_wall, bndval_te_dome, bndval_te_out, & bndval_ti_core, bndval_ti_wall, bndval_ti_dome, bndval_ti_out, & bndval_upar_core, bndval_upar_wall, bndval_upar_dome, bndval_upar_out, & bndval_qe_core, bndval_qe_wall, bndval_qe_dome, bndval_qe_out, & bndval_qi_core, bndval_qi_wall, bndval_qi_dome, bndval_qi_out implicit none type, private :: bndgen_t !! Datatype for defining boundary type and value integer, allocatable, dimension(:) :: types !! Type of boundary condition real(GP), allocatable, dimension(:) :: vals !! Boundary values integer, allocatable, dimension(:) :: solver_types !! Type of boundary condition passed to the elliptic solvers (which accept only DIRICHLET_ZERO or NEUMANN) end type type, public :: boundaries_braginskii_t !! Datatype for boundary types and values of braginskii model type(bndgen_t), public :: ne !! For density type(bndgen_t), public :: pot !! For electrostatic potential type(bndgen_t), public :: vort !! For generalised vorticity type(bndgen_t), public :: apar !! For parallel electromagnetic potential type(bndgen_t), public :: jpar !! For parallel current type(bndgen_t), public :: upar !! For generalised vorticity type(bndgen_t), public :: te !! For electron temperature type(bndgen_t), public :: ti !! For ion temperature type(bndgen_t), public :: qe !! For parallel electron heat flux type(bndgen_t), public :: qi !! For parallel ion heat flux contains procedure, public :: init => init_boundaries_braginskii final :: destructor end type contains subroutine init_boundaries_braginskii(self, mesh_cano, mesh_stag) !! Sets boundary types and values according to parameters class(boundaries_braginskii_t), intent(inout) :: self !! Instance of the type type(mesh_cart_t), intent(in) :: mesh_cano !! Mesh (canonical) type(mesh_cart_t), intent(in) :: mesh_stag !! Mesh (staggered) integer :: nbnds_cano, nbnds_stag, kb, l, district nbnds_cano = mesh_cano%get_n_points_boundary() nbnds_stag = mesh_stag%get_n_points_boundary() ! Boundary descriptors and values for quantities defined on canonical mesh allocate(self%ne%types(nbnds_cano)) allocate(self%ne%vals(nbnds_cano)) allocate(self%pot%types(nbnds_cano)) allocate(self%pot%solver_types(nbnds_cano)) allocate(self%pot%vals(nbnds_cano)) allocate(self%vort%types(nbnds_cano)) allocate(self%vort%vals(nbnds_cano)) allocate(self%te%types(nbnds_cano)) allocate(self%te%solver_types(nbnds_cano)) allocate(self%te%vals(nbnds_cano)) allocate(self%ti%types(nbnds_cano)) allocate(self%ti%solver_types(nbnds_cano)) allocate(self%ti%vals(nbnds_cano)) !$omp parallel default(none) & !$omp private(kb, l, district) & !$omp shared(self, mesh_cano, & !$omp bnddescr_ne_core, bnddescr_ne_wall, bnddescr_ne_dome, bnddescr_ne_out, & !$omp bnddescr_pot_core, bnddescr_pot_wall, bnddescr_pot_dome, bnddescr_pot_out, & !$omp bnddescr_vort_core, bnddescr_vort_wall, bnddescr_vort_dome, bnddescr_vort_out, & !$omp bnddescr_te_core, bnddescr_te_wall, bnddescr_te_dome, bnddescr_te_out, & !$omp bnddescr_ti_core, bnddescr_ti_wall, bnddescr_ti_dome, bnddescr_ti_out, & !$omp bndval_ne_core, bndval_ne_wall, bndval_ne_dome, bndval_ne_out, & !$omp bndval_pot_core, bndval_pot_wall, bndval_pot_dome, bndval_pot_out, & !$omp bndval_vort_core, bndval_vort_wall, bndval_vort_dome, bndval_vort_out, & !$omp bndval_te_core, bndval_te_wall, bndval_te_dome, bndval_te_out, & !$omp bndval_ti_core, bndval_ti_wall, bndval_ti_dome, bndval_ti_out) !$omp do do kb = 1, mesh_cano%get_n_points_boundary() l = mesh_cano%boundary_indices(kb) district = mesh_cano%get_district(l) select case(district) case(DISTRICT_CORE) ! Speciality for zonal Neumann (treated as hybrid of homogeneous Neumann and Dirichlet) if (bnddescr_ne_core == BND_BRAGTYPE_ZONAL_NEUMANN) then self%ne%types(kb) = BND_TYPE_NEUMANN else self%ne%types(kb) = bnddescr_ne_core endif self%ne%vals(kb) = log_bnd(bnddescr_ne_core, bndval_ne_core) if (bnddescr_pot_core == BND_BRAGTYPE_ZONAL_NEUMANN) then ! Specialty for zonal Neumann boundary condition here ! Set boundary to Dirichlet zero, as used in apply_helmholtz from (second) zonal solve ! Shall be improved with new multigrid solver self%pot%types(kb) = BND_TYPE_DIRICHLET_ZERO self%pot%solver_types(kb) = BND_TYPE_DIRICHLET_ZERO else self%pot%types(kb) = bnddescr_pot_core select case (bnddescr_pot_core) case (BND_TYPE_DIRICHLET_ZERO, BND_TYPE_NEUMANN) self%pot%solver_types(kb) = bnddescr_pot_core case default ! Default is to set all unhandled solver types to a dirichlet condition self%pot%solver_types(kb) = BND_TYPE_DIRICHLET_ZERO end select endif self%pot%vals(kb) = bndval_pot_core self%vort%types(kb) = bnddescr_vort_core self%vort%vals(kb) = bndval_vort_core ! Speciality for zonal Neumann (treated as hybrid of homogeneous Neumann and Dirichlet) if (bnddescr_te_core == BND_BRAGTYPE_ZONAL_NEUMANN) then self%te%types(kb) = BND_TYPE_NEUMANN self%te%solver_types(kb) = BND_TYPE_DIRICHLET_ZERO else self%te%types(kb) = bnddescr_te_core endif self%te%vals(kb) = log_bnd(bnddescr_te_core, bndval_te_core) ! Speciality for zonal Neumann (treated as hybrid of homogeneous Neumann and Dirichlet) if (bnddescr_ti_core == BND_BRAGTYPE_ZONAL_NEUMANN) then self%ti%types(kb) = BND_TYPE_NEUMANN self%ti%solver_types(kb) = BND_TYPE_DIRICHLET_ZERO else self%ti%types(kb) = bnddescr_ti_core endif self%ti%vals(kb) = log_bnd(bnddescr_ti_core, bndval_ti_core) case(DISTRICT_WALL) self%ne%types(kb) = bnddescr_ne_wall self%ne%vals(kb) = log_bnd(bnddescr_ne_wall, bndval_ne_wall) self%pot%types(kb) = bnddescr_pot_wall select case (bnddescr_pot_wall) case (BND_TYPE_DIRICHLET_ZERO, BND_TYPE_NEUMANN) self%pot%solver_types(kb) = bnddescr_pot_wall case default ! Default is to set all unhandled solver types to a dirichlet condition self%pot%solver_types(kb) = BND_TYPE_DIRICHLET_ZERO end select self%pot%vals(kb) = bndval_pot_wall self%vort%types(kb) = bnddescr_vort_wall self%vort%vals(kb) = bndval_vort_wall self%te%types(kb) = bnddescr_te_wall self%te%solver_types(kb) = bnddescr_te_wall self%te%vals(kb) = log_bnd(bnddescr_te_wall, bndval_te_wall) self%ti%types(kb) = bnddescr_ti_wall self%ti%solver_types(kb) = bnddescr_ti_wall self%ti%vals(kb) = log_bnd(bnddescr_ti_wall, bndval_ti_wall) case(DISTRICT_DOME) self%ne%types(kb) = bnddescr_ne_dome self%ne%vals(kb) = log_bnd(bnddescr_ne_dome, bndval_ne_dome) self%pot%types(kb) = bnddescr_pot_dome select case (bnddescr_pot_dome) case (BND_TYPE_DIRICHLET_ZERO, BND_TYPE_NEUMANN) self%pot%solver_types(kb) = bnddescr_pot_dome case default ! Default is to set all unhandled solver types to a dirichlet condition self%pot%solver_types(kb) = BND_TYPE_DIRICHLET_ZERO end select self%pot%vals(kb) = bndval_pot_dome self%vort%types(kb) = bnddescr_vort_dome self%vort%vals(kb) = bndval_vort_dome self%te%types(kb) = bnddescr_te_dome self%te%solver_types(kb) = bnddescr_te_dome self%te%vals(kb) = log_bnd(bnddescr_te_dome, bndval_te_dome) self%ti%types(kb) = bnddescr_ti_dome self%ti%solver_types(kb) = bnddescr_ti_dome self%ti%vals(kb) = log_bnd(bnddescr_ti_dome, bndval_ti_dome) case(DISTRICT_OUT) self%ne%types(kb) = bnddescr_ne_out self%ne%vals(kb) = log_bnd(bnddescr_ne_out, bndval_ne_out) self%pot%types(kb) = bnddescr_pot_out select case (bnddescr_pot_out) case (BND_TYPE_DIRICHLET_ZERO, BND_TYPE_NEUMANN) self%pot%solver_types(kb) = bnddescr_pot_out case default ! Default is to set all unhandled solver types to a dirichlet condition self%pot%solver_types(kb) = BND_TYPE_DIRICHLET_ZERO end select self%pot%vals(kb) = bndval_pot_out self%vort%types(kb) = bnddescr_vort_out self%vort%vals(kb) = bndval_vort_out self%te%types(kb) = bnddescr_te_out self%te%solver_types(kb) = bnddescr_te_out self%te%vals(kb) = log_bnd(bnddescr_te_out, bndval_te_out) self%ti%types(kb) = bnddescr_ti_out self%ti%solver_types(kb) = bnddescr_ti_out self%ti%vals(kb) = log_bnd(bnddescr_ti_out, bndval_ti_out) case default !$omp critical call handle_error('District not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('District: ',[district])) !$omp end critical end select enddo !$omp end do !$omp end parallel ! Boundary descriptors and values for quantities defined on staggered mesh allocate(self%apar%types(nbnds_stag)) allocate(self%apar%vals(nbnds_stag)) allocate(self%jpar%types(nbnds_stag)) allocate(self%jpar%vals(nbnds_stag)) allocate(self%upar%types(nbnds_stag)) allocate(self%upar%vals(nbnds_stag)) allocate(self%qe%types(nbnds_stag)) allocate(self%qe%vals(nbnds_stag)) allocate(self%qi%types(nbnds_stag)) allocate(self%qi%vals(nbnds_stag)) !$omp parallel default(none) & !$omp private(kb, l, district) & !$omp shared(self, mesh_stag, & !$omp bnddescr_apar_core, bnddescr_apar_wall, bnddescr_apar_dome, bnddescr_apar_out, & !$omp bnddescr_jpar_core, bnddescr_jpar_wall, bnddescr_jpar_dome, bnddescr_jpar_out, & !$omp bnddescr_upar_core, bnddescr_upar_wall, bnddescr_upar_dome, bnddescr_upar_out, & !$omp bnddescr_qe_core, bnddescr_qe_wall, bnddescr_qe_dome, bnddescr_qe_out, & !$omp bnddescr_qi_core, bnddescr_qi_wall, bnddescr_qi_dome, bnddescr_qi_out, & !$omp bndval_apar_core, bndval_apar_wall, bndval_apar_dome, bndval_apar_out, & !$omp bndval_jpar_core, bndval_jpar_wall, bndval_jpar_dome, bndval_jpar_out, & !$omp bndval_upar_core, bndval_upar_wall, bndval_upar_dome, bndval_upar_out, & !$omp bndval_qe_core, bndval_qe_wall, bndval_qe_dome, bndval_qe_out, & !$omp bndval_qi_core, bndval_qi_wall, bndval_qi_dome, bndval_qi_out) !$omp do do kb = 1, mesh_stag%get_n_points_boundary() l = mesh_stag%boundary_indices(kb) district = mesh_stag%get_district(l) select case(district) case(DISTRICT_CORE) self%apar%types(kb) = bnddescr_apar_core self%apar%vals(kb) = bndval_apar_core self%jpar%types(kb) = bnddescr_jpar_core self%jpar%vals(kb) = bndval_jpar_core self%upar%types(kb) = bnddescr_upar_core self%upar%vals(kb) = bndval_upar_core self%qe%types(kb) = bnddescr_qe_core self%qe%vals(kb) = bndval_qe_core self%qi%types(kb) = bnddescr_qi_core self%qi%vals(kb) = bndval_qi_core case(DISTRICT_WALL) self%apar%types(kb) = bnddescr_apar_wall self%apar%vals(kb) = bndval_apar_wall self%jpar%types(kb) = bnddescr_jpar_wall self%jpar%vals(kb) = bndval_jpar_wall self%upar%types(kb) = bnddescr_upar_wall self%upar%vals(kb) = bndval_upar_wall self%qe%types(kb) = bnddescr_qe_wall self%qe%vals(kb) = bndval_qe_wall self%qi%types(kb) = bnddescr_qi_wall self%qi%vals(kb) = bndval_qi_wall case(DISTRICT_DOME) self%apar%types(kb) = bnddescr_apar_dome self%apar%vals(kb) = bndval_apar_dome self%jpar%types(kb) = bnddescr_jpar_dome self%jpar%vals(kb) = bndval_jpar_dome self%upar%types(kb) = bnddescr_upar_dome self%upar%vals(kb) = bndval_upar_dome self%qe%types(kb) = bnddescr_qe_dome self%qe%vals(kb) = bndval_qe_dome self%qi%types(kb) = bnddescr_qi_dome self%qi%vals(kb) = bndval_qi_dome case(DISTRICT_OUT) self%apar%types(kb) = bnddescr_apar_out self%apar%vals(kb) = bndval_apar_out self%jpar%types(kb) = bnddescr_jpar_out self%jpar%vals(kb) = bndval_jpar_out self%upar%types(kb) = bnddescr_upar_out self%upar%vals(kb) = bndval_upar_out self%qe%types(kb) = bnddescr_qe_out self%qe%vals(kb) = bndval_qe_out self%qi%types(kb) = bnddescr_qi_out self%qi%vals(kb) = bndval_qi_out case default !$omp critical call handle_error('District not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('District: ',[district])) !$omp end critical end select enddo !$omp end do !$omp end parallel end subroutine real(GP) function log_bnd(bnddescr, bndval) !! Sets bnd value to the logarithm of the prescribed value if the !! descrpitor is "DIRICHLET". Necessary for ne, te and ti integer, intent(in) :: bnddescr !! Type of perpendicular boundary condition real(GP), intent(in) :: bndval !! Value for the perpendicular boundary condition if (bnddescr == BND_TYPE_DIRICHLET_ZERO) then log_bnd = log(bndval) else log_bnd = bndval endif end function subroutine zonal_neumann_core(comm_handler, mesh, equi, polars, u, v) !! Sets Zonal Neumann boundary conditions on u (or v), u_b = <u_b_N>, !! after local Neumann b.c. d/dr u_b_N = 0 had been applied to u. !! u is flux-surface averaged and the core value is used as a Dirichlet b.c.. !! For data handling reasons, in- and output can be separated with variable v = <u>. type(comm_handler_t), intent(in) :: comm_handler !! Communicators type(mesh_cart_t), intent(in) :: mesh !! Mesh class(equilibrium_t), intent(in) :: equi !! Equilibrium type(polars_t), intent(in) :: polars !! Polar grid real(GP), dimension(mesh%get_n_points()), intent(inout) :: u !! Variable for which boundary conditions are computed and set real(GP), dimension(mesh%get_n_points()), intent(inout), optional :: v !! Variable for which boundary conditions are set instead of u integer, parameter :: intorder=3 integer :: kb, l, nrho, district real(GP) :: rho_core_rel, drho, u_core real(GP), dimension(polars%grid%get_nrho()) :: u_zon ! Perform zonal average and interpolate to core boundary call surface_average(comm_handler%get_comm_planes(), mesh, polars, u, u_zon) nrho = polars%grid%get_nrho() drho = polars%grid%get_drho() rho_core_rel = equi%rhomin - polars%grid%get_rhopol_min() u_core = interpol1d(intorder, nrho, drho, u_zon, rho_core_rel) !$omp parallel default(none) & !$omp private(kb, l, district) & !$omp shared(mesh, u, u_core, v) !$omp do do kb = 1, mesh%get_n_points_boundary() ! Set boundaries to zonal average at core l = mesh%boundary_indices(kb) district = mesh%get_district(l) if (district == DISTRICT_CORE) then if (present(v)) then v(l) = u_core else u(l) = u_core endif endif enddo !$omp end do !$omp end parallel end subroutine subroutine destructor(self) !! Frees memory associated with boundaries_braginskii type(boundaries_braginskii_t), intent(inout) :: self !! Instance of the type if (allocated(self%ne%types)) deallocate(self%ne%types) if (allocated(self%ne%vals)) deallocate(self%ne%vals) if (allocated(self%pot%types)) deallocate(self%pot%types) if (allocated(self%pot%solver_types)) deallocate(self%pot%solver_types) if (allocated(self%pot%vals)) deallocate(self%pot%vals) if (allocated(self%vort%types)) deallocate(self%vort%types) if (allocated(self%vort%vals)) deallocate(self%vort%vals) if (allocated(self%apar%types)) deallocate(self%apar%types) if (allocated(self%apar%vals)) deallocate(self%apar%vals) if (allocated(self%jpar%types)) deallocate(self%jpar%types) if (allocated(self%jpar%vals)) deallocate(self%jpar%vals) if (allocated(self%upar%types)) deallocate(self%upar%types) if (allocated(self%upar%vals)) deallocate(self%upar%vals) if (allocated(self%te%types)) deallocate(self%te%types) if (allocated(self%te%solver_types)) deallocate(self%te%solver_types) if (allocated(self%te%vals)) deallocate(self%te%vals) if (allocated(self%ti%types)) deallocate(self%ti%types) if (allocated(self%ti%solver_types)) deallocate(self%ti%solver_types) if (allocated(self%ti%vals)) deallocate(self%ti%vals) if (allocated(self%qe%types)) deallocate(self%qe%types) if (allocated(self%qe%vals)) deallocate(self%qe%vals) if (allocated(self%qi%types)) deallocate(self%qi%types) if (allocated(self%qi%vals)) deallocate(self%qi%vals) end subroutine end module