module immersed_rho_m !! immersed boundary module, bsaed on flux surface label rho use precision_m, only : FP, FP_LARGEST use screen_io_m, only : get_stdout use error_handling_m, only : handle_error use status_codes_m, only : PARALLAX_ERR_PARAMETERS use immersed_m, only : immersed_t use equilibrium_m, only : equilibrium_t use mesh_cart_m, only : mesh_cart_t use elementary_functions_m, only : step_hermite implicit none type, public, extends(immersed_t) :: immersed_rho_t !! Immersed boundary type based on rho contains procedure, public :: init procedure, public :: display procedure, private :: read_params_immersed_rho final :: destructor_immersed_rho_t end type ! Parameters related with rho immersed boundary real(FP), protected :: immersed_rho_inner = -FP_LARGEST !! rho value, until where boundary immersion is applied real(FP), protected :: immersed_width_inner = 0.0_FP !! Transition width at inner rho boundary immersion real(FP), protected :: immersed_rho_outer = FP_LARGEST !! rho value, from where boundary immersion is applied real(FP), protected :: immersed_width_outer = 0.0_FP !! Transition width at outer rho boundary immersion integer, protected :: immersed_step_order = 2 !! Hermite order of transition function namelist / immersed_rho / & immersed_rho_inner, & immersed_width_inner, & immersed_rho_outer, & immersed_width_outer, & immersed_step_order contains subroutine init(self, equi, mesh, filename) class(immersed_rho_t), intent(inout) :: self class(equilibrium_t), intent(inout) :: equi type(mesh_cart_t), intent(in) :: mesh character(len=*), intent(in), optional :: filename integer :: l real(FP) :: x, y, phi, rho if (present(filename)) then call self%read_params_immersed_rho(filename) endif phi = mesh%get_phi() ! Allocate fields allocate(self%charfun(mesh%get_n_points())) allocate(self%dirindfun(mesh%get_n_points())) !$omp parallel default(none) private(l, x, y, rho) & !$omp shared(self, mesh, equi, phi, immersed_rho_outer, & !$omp immersed_width_outer, immersed_rho_inner, & !$omp immersed_width_inner, immersed_step_order) !$omp do do l = 1, mesh%get_n_points() x = mesh%get_x(l) y = mesh%get_y(l) rho = equi%rho(x, y, phi) self%charfun(l) = step_hermite(immersed_rho_outer, rho, & immersed_width_outer, order=immersed_step_order) & + 1.0_FP - step_hermite(immersed_rho_inner, rho, & immersed_width_inner, order=immersed_step_order) ! Set dirindfun to 0 self%dirindfun(l) = 0.0_FP enddo !$omp end do !$omp end parallel ! Build inds and adj_inds call self%build_inds(mesh) call self%build_adj_inds(mesh) end subroutine subroutine read_params_immersed_rho(self, filename) !! Reads parameters for rho based immersed boundary class(immersed_rho_t), intent(inout) :: self !! Instance of type character(len=*), intent(in):: filename !! Filename, where to read parameters from integer :: io_unit, io_error character(len=256) :: io_errmsg open(newunit=io_unit, file=filename, status='old', action='read', & iostat=io_error, iomsg=io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, PARALLAX_ERR_PARAMETERS, & __LINE__, __FILE__) endif read(io_unit, nml=immersed_rho, iostat=io_error, iomsg=io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, PARALLAX_ERR_PARAMETERS, & __LINE__, __FILE__) endif close(io_unit) end subroutine subroutine display(self) class(immersed_rho_t), intent(in) :: self integer :: io_unit, io_error character(len=256) :: io_errmsg io_unit = get_stdout() write(io_unit, nml=immersed_rho, iostat=io_error, iomsg=io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, PARALLAX_ERR_PARAMETERS, & __LINE__, __FILE__) endif end subroutine subroutine destructor_immersed_rho_t(self) !! Destructor for immersed_rho_t type(immersed_rho_t), intent(inout) :: self if (allocated(self%charfun)) deallocate(self%charfun) if (allocated(self%dirindfun)) deallocate(self%dirindfun) if (allocated(self%inds)) deallocate(self%inds) if (allocated(self%adj_inds)) deallocate(self%adj_inds) end subroutine end module