module sources_external_m !! External sources use MPI use comm_handler_m, only : comm_handler_t use precision_grillix_m, only : GP, GP_EPS use constants_m, only : TWO_PI use error_handling_grillix_m, only: handle_error, error_info_t use status_codes_grillix_m, only : GRILLIX_ERR_OTHER, GRILLIX_ERR_NAMELIST use screen_io_m, only : get_stdout use parallelisation_setup_m, only : is_rank_info_writer use variable_m, only : variable_t use equilibrium_storage_m, only : equilibrium_storage_t use perf_m, only : perf_start, perf_stop use parallel_map_m, only : parallel_map_t use mesh_cart_m, only : mesh_cart_t use polars_m, only : polars_t use polar_operators_m, only : surface_average use elementary_functions_m, only : gaussian, step_tanh, step_hermite, gaussian use interpolation_m, only : interpol1d use descriptors_m, only : DISTRICT_CLOSED use technical_constants_m, only : PATHLEN_MAX use diagnostics_helpers_m, only : inner_product use penalisation_m, only : penalisation_t use snapshot_m, only : snapshot_t, char_arr_t use params_neut_model_m, only : gas_puff_on, gas_puff_path ! Parameters use params_brag_model_m, only : & tratio use params_brag_sources_external_select_m, only : & nesrc_select, nesrc_path, tesrc_select, tesrc_path, & tisrc_select, tisrc_path, is_powersource implicit none type, private, abstract :: source_individual_t !! Abstract class for an individual source contains procedure(init_source_individual), deferred :: init procedure(display_source_individual), deferred :: display end type interface subroutine init_source_individual(self, filename) import source_individual_t !!Initialises an individual source class(source_individual_t), intent(inout) :: self !! Instance of type character(len=*), intent(in) :: filename !! Filename to read parameters from end subroutine subroutine display_source_individual(self) import source_individual_t !! Displays parameters of an individual source class(source_individual_t), intent(in) :: self !! Instance of type end subroutine end interface type, public, extends(source_individual_t) :: source_none_t !! Parameters for none source contains procedure :: init => init_source_none procedure :: display => display_source_none end type type, public, extends(source_individual_t) :: source_zonal_t !! Zonal adaptive source !! Damps zonal averaged quantities to prescribed profile character(len=32), private :: srcfun_type !! Profile type procedure(srcfun_interface), private, nopass, pointer :: srcfun => null() !! Pointer to function of source profile real(GP), private :: rho_loc = 0.0_GP !! Location of source real(GP) :: width = 1.0_GP !! Width of source (rho space) real(GP), private :: rate = 0.0_GP !! Source rate (strength) integer, private :: n_profile !! Number of points of prescribed profile real(GP), private :: rho_ref_profile !! Reference point of (equidistant) profile grid (smallest rho) real(GP), private :: drho_profile !! Distance of (equidistant) profile grid !! rho_profile(i) = rho_ref_profile + (i-1) * drho_profile, !! i = 1, n_profile real(GP), private, allocatable, dimension(:) :: vals_profile !! Values of prescribed profile contains procedure :: eval_source_zonal procedure :: init => init_source_zonal procedure :: display => display_source_zonal end type type, public, extends(source_individual_t) :: source_constantrate_t !! Source with constant rate (non adaptive) character(len=32), private :: srcfun_type !! Profile type procedure(srcfun_interface), private, nopass, pointer :: srcfun => null() !! Pointer to function of source profile real(GP), private :: rho_loc !! Location of source real(GP), private :: width !! Width of source (rho space) real(GP), private :: rate !! Source rate (strength) character(len=32), private :: read_from !! Indicates from which parameterfile source was read contains procedure :: eval_source_constantrate procedure :: init => init_source_constantrate procedure :: display => display_source_constantrate end type private srcfun_interface abstract interface real(GP) function srcfun_interface(rho, rho_loc, wrho) !! Interface for source profile function import GP real(GP), intent(in) :: rho real(GP), intent(in) :: rho_loc real(GP), intent(in) :: wrho end function end interface private :: srcfun_gauss private :: srcfun_step_tanh private :: srcfun_step_hermite private :: srcfun_tcv_special type, public :: sources_external_t !! External sources class(source_individual_t), allocatable, private :: par_srcne !! Parameter for external particle source class(source_individual_t), allocatable, private :: par_srcte !! Parameter for external electron temperature source class(source_individual_t), allocatable, private :: par_srcti !! Parameter for external ion temperature source contains procedure, public :: init => init_sources_external procedure, public :: eval => eval_sources_external end type type, public, extends(source_individual_t) :: source_gaussian_t !! Class for gaussian sources, contains set of individual source points logical, public :: source_gaussian_on !! Switch whether or not to use gaussian source module character(len=:), allocatable, private :: source_gaussian_path !! Path to file where source parameters are specified integer, private :: n_sources !! Number of individual sourcepoints real(GP), dimension(:), allocatable, private :: x !! x-coordinates of the center of the source points real(GP), dimension(:), allocatable, private :: y !! y-coordinates of the center of the source points real(GP), dimension(:), allocatable, private :: r_width !! Radial widths of the gaussian sources in the plane real(GP), dimension(:), allocatable, private :: phi !! Toroidal angles at the center of the source points real(GP), dimension(:), allocatable, private :: phi_width !! Toroidal widths of the gaussian sources real(GP), dimension(:), allocatable, private :: rate !! Sourcing rates in units u_0/t_0 real(GP), dimension(:), allocatable, public :: src_coeffs !! Precomputed coefficients used for computing localised sources during timestepping contains procedure :: init => init_source_gaussian procedure :: set_coeffs procedure :: display => display_source_gaussian procedure :: eval_source_gaussian final :: destructor end type contains subroutine init_source_gaussian(self, filename) !! Initialize gaussian source from input file class(source_gaussian_t), intent(inout) :: self !! Instance of class character(len=*), intent(in) :: filename !! Path to input file integer :: i, io_unit, io_error character(len=256) :: io_errmsg character(len=1024) :: line self%source_gaussian_on = gas_puff_on self%source_gaussian_path = trim(adjustl(gas_puff_path)) io_unit = 29 ! Scan parameter file if ( self%source_gaussian_on ) then open(unit = io_unit, file = self%source_gaussian_path, status = 'old', & action = 'read', iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif self%n_sources = 0 ! Read line by line and count non-comment lines do read(io_unit, '(a)', iostat=io_error) line if (io_error /= 0) exit if ((line(1:1) == '!') .or. (line == "") .or. (line(1:1) == " ")) cycle ! Skip comment lines self%n_sources = self%n_sources + 1 enddo rewind(io_unit) ! Restart from top allocate(self%x(self%n_sources)) allocate(self%y(self%n_sources)) allocate(self%r_width(self%n_sources)) allocate(self%phi(self%n_sources)) allocate(self%phi_width(self%n_sources)) allocate(self%rate(self%n_sources)) i = 1 do while ( i <= self%n_sources ) read(io_unit, '(a)', iostat=io_error) line if (io_error /= 0) exit if ((line(1:1) == '!') .or. (line == "") .or. (line(1:1) == " ")) cycle ! Skip comment lines ! Read source point data read(line, *, iostat = io_error, iomsg = io_errmsg) & self%x(i), self%y(i), & self%r_width(i), self%phi(i), & self%phi_width(i), self%rate(i) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif i = i + 1 enddo close(io_unit) endif end subroutine subroutine set_coeffs(self, comm_handler, mesh_cano, map, penalisation_cano) !! Compute the coefficient that matches the total input rate class(source_gaussian_t), intent(inout) :: self !! Instance of class type(comm_handler_t), intent(in) :: comm_handler !! Communicators type(mesh_cart_t), intent(inout) :: mesh_cano !! Mesh (canonical) within poloidal plane type(parallel_map_t), intent(in) :: map !! Parallel map type(penalisation_t), intent(in) :: penalisation_cano !! Penalisation (canonical) real(GP), dimension(mesh_cano%get_n_points()) :: source_on_mesh real(GP) :: total_source integer :: i ! Initialize with 1.0, as some values are required for evaluating gaussian sources allocate(self%src_coeffs(self%n_sources), source=1.0_GP) do i = 1, self%n_sources call self%eval_source_gaussian(mesh_cano, map, source_on_mesh, i_source=i) total_source = inner_product(comm_handler, mesh_cano, map, & penalisation_cano, is_stag=.false., u_in=source_on_mesh, & use_abs=.false., & use_vol_avg=.false., & use_full_domain=.true., & sum_over_planes=.true.) self%src_coeffs(i) = self%rate(i) / total_source enddo end subroutine subroutine eval_source_gaussian(self, mesh, map, source_on_mesh, i_source) !! Compute neutrals density source from contributions of gaussian centers class(source_gaussian_t), intent(inout) :: self !! Instance of class type(mesh_cart_t), intent(inout) :: mesh !! Mesh (canonical or staggered) within poloidal plane type(parallel_map_t), intent(in) :: map !! Parallel map real(GP), dimension(mesh%get_n_points()), intent(inout) :: source_on_mesh !! Results of the source on the each mesh point integer, optional :: i_source !! Index of gaussian source to evaluate. If not provided, evaluates sum of !! all gaussian sources real(GP) :: co_toroidal, delta_phi integer :: i, l, i_start, i_end if ( present(i_source) ) then i_start = i_source i_end = i_source else i_start = 1 i_end = self%n_sources endif source_on_mesh = 0.0_GP do i = i_start, i_end ! Compute the toroidal distance between source point and current plane delta_phi = min(TWO_PI - abs(self%phi(i) - mesh%get_phi()), abs(self%phi(i) - mesh%get_phi())) if ( self%phi_width(i) > GP_EPS) then co_toroidal = gaussian(0.0_GP, self%phi_width(i), delta_phi) else if ( delta_phi < map%get_dphi() / 2.0_GP ) then co_toroidal = 1.0_GP else co_toroidal = 0.0_GP endif endif ! Compute the inplane distribution !$omp parallel default(none) private(l) & !$omp shared(mesh, self, source_on_mesh, co_toroidal, i) !$omp do do l = 1, mesh%get_n_points() source_on_mesh(l) = source_on_mesh(l) & + co_toroidal * gaussian(self%x(i), self%r_width(i), mesh%get_x(l)) & * gaussian(self%y(i), self%r_width(i), mesh%get_y(l)) & * self%src_coeffs(i) enddo !$omp end do !$omp end parallel enddo end subroutine subroutine display_source_gaussian(self) !! Display parameters related to the gas puff module class(source_gaussian_t), intent(in) :: self !! Instance of class integer :: i, io_error character(len=256) :: io_errmsg if (.not.is_rank_info_writer) return write(get_stdout(),2103)'gaussian src #', 'x', 'y', 'r_width', 'phi', 'phi_width', 'rate', 'src_coeffs' 2103 FORMAT(A15, A10, A10, A10, A10, A10, A12, A12) do i = 1, self%n_sources write(get_stdout(), 2104, iostat = io_error, iomsg = io_errmsg) & i, self%x(i), self%y(i), & self%r_width(i), self%phi(i), & self%phi_width(i), self%rate(i), & self%src_coeffs(i) 2104 FORMAT(I15, F10.3, F10.3, F10.3, F10.3, F10.3, ES12.3E3, ES12.3E3) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif enddo end subroutine subroutine destructor(self) !! Destructor type(source_gaussian_t), intent(inout) :: self !! Instance of the type self%n_sources = 0 if (allocated(self%x)) deallocate(self%x) if (allocated(self%y)) deallocate(self%y) if (allocated(self%r_width)) deallocate(self%r_width) if (allocated(self%phi)) deallocate(self%phi) if (allocated(self%phi_width)) deallocate(self%phi_width) if (allocated(self%rate)) deallocate(self%rate) if (allocated(self%src_coeffs)) deallocate(self%src_coeffs) end subroutine subroutine init_sources_external(self) !! Initialises (reads parameters) sources_external class(sources_external_t), intent(inout) :: self !! Instance of class type(source_zonal_t), allocatable :: source_zonal type(source_constantrate_t), allocatable :: source_constantrate type(source_none_t), allocatable :: source_none ! Allocate type of source parameters select case(nesrc_select) case('source_zonal') allocate(source_zonal) call move_alloc(source_zonal, self%par_srcne) case('source_constantrate') allocate(source_constantrate) call move_alloc(source_constantrate, self%par_srcne) case default allocate(source_none) call move_alloc(source_none, self%par_srcne) end select select case(tesrc_select) case('source_zonal') allocate(source_zonal) call move_alloc(source_zonal, self%par_srcte) case('source_constantrate') allocate(source_constantrate) call move_alloc(source_constantrate, self%par_srcte) case default allocate(source_none) call move_alloc(source_none, self%par_srcte) end select select case(tisrc_select) case('source_zonal') allocate(source_zonal) call move_alloc(source_zonal, self%par_srcti) case('source_constantrate') allocate(source_constantrate) call move_alloc(source_constantrate, self%par_srcti) case default allocate(source_none) call move_alloc(source_none, self%par_srcti) end select ! Read parameters call self%par_srcne%init(nesrc_path) call self%par_srcte%init(tesrc_path) call self%par_srcti%init(tisrc_path) ! Display parameters if (is_rank_info_writer) then write(get_stdout(),*)'' write(get_stdout(),*)'Particle source parameters -----------------------------' call self%par_srcne%display() write(get_stdout(),*)'--------------------------------------------------------' write(get_stdout(),*)'' write(get_stdout(),*)'' write(get_stdout(),*)'Electron temperature source parameters -----------------' call self%par_srcte%display() write(get_stdout(),*)'--------------------------------------------------------' write(get_stdout(),*)'' write(get_stdout(),*)'' write(get_stdout(),*)'Ion temperature source parameters ----------------------' call self%par_srcti%display() write(get_stdout(),*)'--------------------------------------------------------' write(get_stdout(),*)'' endif end subroutine subroutine eval_sources_external(self, comm_handler, equi_storage, mesh, polars, tau, & ne, pot, vort, upar, jpar, apar, te, ti, & src_ne, src_te, src_ti ) !! Evaluates external sources type(comm_handler_t), intent(in) :: comm_handler !! Communicators class(sources_external_t), intent(inout) :: self !! Instance of class class(equilibrium_storage_t), intent(in) :: equi_storage !! Equilibrium storage type(mesh_cart_t), intent(in) :: mesh !! Mesh within poloidal plane type(polars_t), intent(in) :: polars !! Polar grid and operators real(GP), intent(in) :: tau !! Time, at output tau + dtau type(variable_t), intent(in) :: ne !! Electron density type(variable_t), intent(in) :: pot !! Electrostatic potential type(variable_t), intent(in) :: vort !! Generalised vorticity type(variable_t), intent(in) :: upar !! Parallel ion velocity type(variable_t), intent(in) :: jpar !! Parallel current type(variable_t), intent(in) :: apar !! Parallel electromagnetic potential type(variable_t), intent(in) :: te !! Electron temperature type(variable_t), intent(in) :: ti !! Ion temperature real(GP), dimension(mesh%get_n_points_inner()), intent(inout) :: src_ne !! Particel source on inner grid points real(GP), dimension(mesh%get_n_points_inner()), intent(inout) :: src_te !! Electron temperature source on inner grid points real(GP), dimension(mesh%get_n_points_inner()), intent(inout) :: src_ti !! Ion temperature source on inner grid points integer :: ki, l real(GP), dimension(polars%grid%get_nrho()) :: wrk_zon real(GP), dimension(mesh%get_n_points_inner()):: src_wrk class(source_individual_t), pointer :: tmp => null() call perf_start('sources_external') ! Particle source select type(tmp => self%par_srcne) type is(source_zonal_t) call surface_average(comm_handler%get_comm_planes(), mesh, polars, ne%vals, wrk_zon) call tmp%eval_source_zonal(equi_storage, mesh, polars, wrk_zon, src_ne) type is(source_constantrate_t) call tmp%eval_source_constantrate(equi_storage, mesh, src_ne) end select ! Electron temperature source !$omp parallel default(none) private(ki) shared(mesh, src_wrk) !$omp do do ki = 1, mesh%get_n_points_inner() src_wrk(ki) = 0.0_GP enddo !$omp end do !$omp end parallel select type(tmp => self%par_srcte) type is(source_zonal_t) call surface_average(comm_handler%get_comm_planes(), mesh, polars, te%vals, wrk_zon) call tmp%eval_source_zonal(equi_storage, mesh, polars, wrk_zon, src_wrk) type is(source_constantrate_t) call tmp%eval_source_constantrate(equi_storage, mesh, src_wrk) end select if (is_powersource) then !$omp parallel default(none) private(ki, l) shared(mesh, src_wrk, src_te, src_ne, ne, te) !$omp do do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) src_te(ki) = src_te(ki) + (2.0_GP/3.0_GP * src_wrk(ki) - te%vals(l) * src_ne(ki)) / ne%vals(l) enddo !$omp end do !$omp end parallel else !$omp parallel default(none) private(ki) shared(mesh, src_wrk, src_te) !$omp do do ki = 1, mesh%get_n_points_inner() src_te(ki) = src_te(ki) + src_wrk(ki) enddo !$omp end do !$omp end parallel endif ! Ion temperature source !$omp parallel default(none) private(ki) shared(mesh, src_wrk) !$omp do do ki = 1, mesh%get_n_points_inner() src_wrk(ki) = 0.0_GP enddo !$omp end do !$omp end parallel select type(tmp => self%par_srcti) type is(source_zonal_t) call surface_average(comm_handler%get_comm_planes(), mesh, polars, ti%vals, wrk_zon) call tmp%eval_source_zonal(equi_storage, mesh, polars, wrk_zon, src_wrk) type is(source_constantrate_t) call tmp%eval_source_constantrate(equi_storage, mesh, src_wrk) end select if (is_powersource) then !$omp parallel default(none) private(ki, l) shared(mesh, src_wrk, src_ti, src_ne, ne, ti, tratio) !$omp do do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) src_ti(ki) = src_ti(ki) + (2.0_GP/3.0_GP * src_wrk(ki) - tratio * ti%vals(l) * src_ne(ki)) / ne%vals(l) enddo !$omp end do !$omp end parallel else !$omp parallel default(none) private(ki) shared(mesh, src_wrk, src_ti) !$omp do do ki = 1, mesh%get_n_points_inner() src_ti(ki) = src_ti(ki) + src_wrk(ki) enddo !$omp end do !$omp end parallel endif call perf_stop('sources_external') end subroutine subroutine init_source_none(self, filename) !! Initialises source=none class(source_none_t), intent(inout) :: self !! Instance of type character(len=*), intent(in) :: filename !! Filename to read parameters from ! Nothing to do end subroutine subroutine display_source_none(self) !! Displays source_none parameters class(source_none_t), intent(in) :: self !! Instance of type if (.not.is_rank_info_writer) return write(get_stdout(),*)'Source = NONE' end subroutine subroutine init_source_zonal(self, filename) !! Initialises zonal source class(source_zonal_t), intent(inout) :: self !! Instance of type character(len=*), intent(in) :: filename !! Filename to read parameters from integer, parameter :: n_profile_max = 256 !! Maximum number of profile points that can be prescribed integer :: n_profile real(GP) :: rho_loc, width, rate, rho_ref_profile, drho_profile real(GP), dimension(n_profile_max) :: vals_profile character(len=32) :: srcfun_type integer :: io_error character(len=256) :: io_errmsg namelist / source_zonal / srcfun_type, rho_loc, width, rate, & n_profile, rho_ref_profile, drho_profile, vals_profile open(unit = 25, file = filename, status = 'old', action = 'read',& iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif ! default values srcfun_type = 'GAUSS' rho_loc = 0.0_GP width = 1.0_GP rate = 0.0_GP n_profile = 1 rho_ref_profile = 0.0_GP drho_profile = 0.0_GP vals_profile = 0.0_GP read(25, nml = source_zonal, iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif self%srcfun_type = srcfun_type self%rho_loc = rho_loc self%width = width self%rate = rate self%n_profile = n_profile self%rho_ref_profile = rho_ref_profile if ( (self%n_profile > 1) .and. (self%n_profile < 4)) then call handle_error('N_profile must be 1 or larger equal 4', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('N_profile: ',[self%n_profile])) endif self%drho_profile = drho_profile allocate(self%vals_profile(self%n_profile)) self%vals_profile(1:self%n_profile) = vals_profile(1:self%n_profile) select case(self%srcfun_type) case('GAUSS') self%srcfun => srcfun_gauss case('STEP_TANH') self%srcfun => srcfun_step_tanh case('STEP_HERMITE') self%srcfun => srcfun_step_hermite case('TCV_SPECIAL') self%srcfun => srcfun_tcv_special case default call handle_error('Srcfun_type not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t(self%srcfun_type)) end select close(25) end subroutine subroutine display_source_zonal(self) !! Displays zonal source parameters class(source_zonal_t), intent(in) :: self !! Instance of type integer :: iprf if (.not.is_rank_info_writer) return write(get_stdout(),*)'Zonal source:' write(get_stdout(),101)self%srcfun_type, self%rho_loc, self%width, self%rate,& self%n_profile, self%rho_ref_profile, self%drho_profile 101 FORMAT(3X,'srcfun_type = ',A32 /, & 3X,'rho_loc = ',ES14.6E3 /, & 3X,'width = ',ES14.6E3 /, & 3X,'rate = ',ES14.6E3 /, & 3X,'n_profile = ',I4 /, & 3X,'rho_ref_profile = ',ES14.6E3 /, & 3X,'drho_profile = ',ES14.6E3 ) do iprf = 1, self%n_profile write(get_stdout(),102)self%rho_ref_profile + (iprf - 1) & * self%drho_profile, self%vals_profile(iprf) 102 FORMAT(4X,'rho_profile = ',X,ES14.6E3,X,', val_profile = ',X,ES14.6E3) enddo end subroutine subroutine eval_source_zonal(self, equi_storage, mesh, polars, uzon, src_u) !! Evaluates zonal source class(source_zonal_t), intent(in) :: self !! Instance of type class(equilibrium_storage_t), intent(in) :: equi_storage !! Equilibrium storage type(mesh_cart_t), intent(in) :: mesh !! Mesh type(polars_t), intent(in) :: polars !! Polars real(GP), dimension(polars%grid%get_nrho()), intent(in) :: uzon !! Zonal values of quantity u that is sourced real(GP), dimension(mesh%get_n_points_inner()), intent(inout) :: src_u !! Source values on inner mesh points integer, parameter :: intorder=3 integer :: ki, l, district real(GP) :: x, y, rho, rho_rel, fac, val_prof_current, val_prof_target !$omp parallel default(none) & !$omp shared(self, equi_storage, mesh, polars, uzon, src_u) & !$omp private(ki, l, x, y, rho, rho_rel, fac, & !$omp val_prof_current, val_prof_target, district) !$omp do do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) x = mesh%get_x(l) y = mesh%get_y(l) rho = equi_storage%rho(l) district = equi_storage%district(l) if (district == DISTRICT_CLOSED) then ! Interpolate current profile (zonal average) to local position rho_rel = rho - polars%grid%get_rhopol_min() val_prof_current = interpol1d(intorder, polars%grid%get_nrho(), & polars%grid%get_drho(), & uzon, & rho_rel) ! Get target profile at local position if (self%n_profile == 1) then val_prof_target = self%vals_profile(1) else rho_rel = rho - self%rho_ref_profile val_prof_target = interpol1d(intorder, self%n_profile, & self%drho_profile, & self%vals_profile, & rho_rel) endif fac = self%rate*self%srcfun(rho, self%rho_loc, self%width) src_u(ki) = src_u(ki) & + fac * (val_prof_target - val_prof_current ) endif enddo !$omp end do !$omp end parallel end subroutine subroutine init_source_constantrate(self, filename) !! Initialises constantrate source class(source_constantrate_t), intent(inout) :: self !! Instance of type character(len=*), intent(in) :: filename !! Filename to read parameters from integer :: io_error character(len=256) :: io_errmsg real(GP) :: rho_loc, width, rate character(len=32) :: srcfun_type namelist / source_constantrate / srcfun_type, rho_loc, width, rate open(unit = 25, file = filename, status = 'old', action = 'read',& iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif self%read_from = filename ! default values srcfun_type = 'GAUSS' rho_loc = 0.0_GP width = 1.0_GP rate = 0.0_GP read(25, nml = source_constantrate, iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif self%srcfun_type = srcfun_type self%rho_loc = rho_loc self%width = width self%rate = rate select case(self%srcfun_type) case('GAUSS') self%srcfun => srcfun_gauss case('STEP_TANH') self%srcfun => srcfun_step_tanh case('STEP_HERMITE') self%srcfun => srcfun_step_hermite case('TCV_SPECIAL') self%srcfun => srcfun_tcv_special case default call handle_error('Srcfun_type not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t(self%srcfun_type)) end select close(25) end subroutine subroutine display_source_constantrate(self) !! Displays constantrate source parameters class(source_constantrate_t), intent(in) :: self !! Instance of type if (.not.is_rank_info_writer) return write(get_stdout(),*)'Constantrate source:' write(get_stdout(),102)self%srcfun_type, self%rho_loc, self%width, self%rate 102 FORMAT(3X,'srcfun_type = 'A32 /, & 3X,'rho_loc = ',ES14.6E3 /, & 3X,'width = ',ES14.6E3 /, & 3X,'rate = ',ES14.6E3 ) end subroutine subroutine eval_source_constantrate(self, equi_storage, mesh, src_u) !! Evaluates constantrate source class(source_constantrate_t), intent(in) :: self !! Instance of type class(equilibrium_storage_t), intent(in) :: equi_storage !! Equilibrium storage type(mesh_cart_t), intent(in) :: mesh !! Mesh real(GP), dimension(mesh%get_n_points_inner()), intent(inout) :: src_u !! Source values on inner mesh points integer :: ki, l, district real(GP) :: x, y, rho !$OMP PARALLEL PRIVATE(ki, l, x, y, rho, district) !$OMP DO do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) x = mesh%get_x(l) y = mesh%get_y(l) rho = equi_storage%rho(l) district = equi_storage%district(l) if (district == DISTRICT_CLOSED) then src_u(ki) = src_u(ki) + self%rate * self%srcfun(rho, self%rho_loc, self%width) endif enddo !$OMP END DO !$OMP END PARALLEL end subroutine real(GP) function srcfun_gauss(rho, rho_loc, wrho) !! Source function with Gauss profile real(GP), intent(in) :: rho !! Rho real(GP), intent(in) :: rho_loc !! Center of Gaussian real(GP), intent(in) :: wrho !! Width of Gaussian srcfun_gauss = gaussian(x0 = rho_loc, wx = wrho, x = rho) end function real(GP) function srcfun_step_tanh(rho, rho_loc, wrho) !! Source function with decaying tanh profile real(GP), intent(in) :: rho !! Rho real(GP), intent(in) :: rho_loc !! Center of step real(GP), intent(in) :: wrho !! Width of step srcfun_step_tanh = 1.0_GP - step_tanh(x0 = rho_loc, x = rho, wx = wrho) end function real(GP) function srcfun_step_hermite(rho, rho_loc, wrho) !! Source function with decaying hermite profile real(GP), intent(in) :: rho !! Rho real(GP), intent(in) :: rho_loc !! Center of step real(GP), intent(in) :: wrho !! Width of step srcfun_step_hermite = 1.0_GP - step_hermite(x0 = rho_loc, x = rho, wx = wrho) end function real(GP) function srcfun_tcv_special(rho, rho_loc, wrho) !! Source function for density sourcing of TCV real(GP), intent(in) :: rho !! Rho real(GP), intent(in) :: rho_loc !! Center of step real(GP), intent(in) :: wrho !! Width of step srcfun_tcv_special = exp(-((rho**2 - rho_loc**2)/wrho)**2) end function end module