module immersed_trace_m !! Immersed boundary module, based on parallel distance to targets (tracing) use precision_m, only : FP use screen_io_m, only : get_stdout use constants_m, only : PI use error_handling_m, only : handle_error use status_codes_m, only : PARALLAX_ERR_PARAMETERS, & PARALLAX_ERR_IMMERSED, PARALLAX_WRN_IMMERSED use immersed_m, only : immersed_t use equilibrium_m, only : equilibrium_t use mesh_cart_m, only : mesh_cart_t use connection_length_m, only : distance_to_boundary use elementary_functions_m, only : step_hermite use root_finding_m, only : func_1d_t, find_zero implicit none type, public, extends(immersed_t) :: immersed_trace_t !! Immersed boundary type based on field line tracing contains procedure, public :: init procedure, public :: display procedure, private :: read_params_immersed_trace final :: destructor_immersed_trace_t end type ! Parameters related with trace boundary immersion real(FP), protected :: charfun_at_target = 0.5_FP !! Value of characteristic immersion function at target real(FP), protected :: dphi_max = PI !! Maximum distance of tracing integer, protected :: immersed_step_order = 2 !! Hermite order of transition function real(FP), protected :: immersed_trace_width_charfun = 0.0_FP !! Parallel transition width of characteristic immersion function real(FP), protected :: immersed_trace_width_dirindfun = 0.0_FP !! Parallel transition width length of direction indicator function real(FP), protected :: trace_max_step_size = PI / 180.0_FP !! Maximum step size of trace algorithm type, private, extends(func_1d_t) :: aux_sigmoid_func_t !! Auxiliary function for setting charfun at target contains procedure :: func => aux_sigmoid_func end type namelist / immersed_trace / & charfun_at_target, & dphi_max, & immersed_step_order, & immersed_trace_width_charfun, & immersed_trace_width_dirindfun, & trace_max_step_size contains subroutine init(self, equi, mesh, filename) class(immersed_trace_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) :: phi, x, y, arclength_fwd, arclength_bwd, & dphi_fwd, dphi_bwd, charfun_fwd, charfun_bwd, deltaphi_rel, & dphi_shift logical :: in_vessel, target_reached_fwd, target_reached_bwd if (present(filename)) then call self%read_params_immersed_trace(filename) endif dphi_shift = get_dphi_shift() allocate(self%charfun(mesh%get_n_points())) allocate(self%dirindfun(mesh%get_n_points())) phi = mesh%get_phi() !$omp parallel default(none) & !$omp private(l, x, y, in_vessel, arclength_fwd, arclength_bwd, & !$omp target_reached_fwd, target_reached_bwd, & !$omp dphi_fwd, dphi_bwd, charfun_fwd, charfun_bwd, & !$omp deltaphi_rel) & !$omp shared(self, mesh, equi, phi, dphi_max, dphi_shift, & !$omp immersed_trace_width_charfun, trace_max_step_size, & !$omp immersed_trace_width_dirindfun, immersed_step_order) !$omp do do l = 1, mesh%get_n_points() x = mesh%get_x(l) y = mesh%get_y(l) in_vessel = equi%in_vessel(x, y, phi) call distance_to_boundary(equi, x, y, phi, dphi_max, & arclength_fwd, arclength_bwd, & target_reached_fwd, target_reached_bwd, trace_max_step_size, & dphi_fwd, dphi_bwd) ! Closed field line region or region inside wall if ((.not. target_reached_fwd) .and. & (.not. target_reached_bwd)) then if (in_vessel) then self%charfun(l) = 0.0_FP self%dirindfun(l) = 0.0_FP else self%charfun(l) = 1.0_FP self%dirindfun(l) = 0.0_FP call handle_error('trace_max_step_size for immersed & boundary probably not sufficient, check charfun and & dirindfun and rerun with increased & trace_max_step_size', & PARALLAX_WRN_IMMERSED, __LINE__, __FILE__) endif cycle endif ! Characteristic function charfun_fwd = sigmoid(dphi_shift, dphi_fwd, & immersed_trace_width_charfun, immersed_step_order) charfun_bwd = sigmoid(dphi_shift, dphi_bwd, & immersed_trace_width_charfun, immersed_step_order) if (abs(dphi_bwd) > abs(dphi_fwd)) then self%charfun(l) = charfun_fwd else self%charfun(l) = charfun_bwd endif ! Direction indicator function deltaphi_rel = abs(dphi_bwd) - abs(dphi_fwd) self%dirindfun(l) = 1.0_FP - & 2.0_FP * sigmoid(0.0_FP, deltaphi_rel, & immersed_trace_width_dirindfun, immersed_step_order) 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_trace(self, filename) !! Reads parameters for trace based immersion class(immersed_trace_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_trace, 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_trace_t), intent(in) :: self integer :: io_unit, io_error character(len=256) :: io_errmsg io_unit = get_stdout() write(io_unit, nml=immersed_trace, iostat=io_error, iomsg=io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, PARALLAX_ERR_PARAMETERS, & __LINE__, __FILE__) endif end subroutine real(FP) function sigmoid(x0, x, wx, order) !! Sigmoid function going from y=0 to y=1 between x=x0-wx/2 to x=x0+wx/2 real(FP), intent(in) :: x0 !! Location of center of sigmoid real(FP), intent(in) :: x !! Abscissa real(FP), intent(in) :: wx !! Transition width integer, intent(in) :: order !! Order of transition sigmoid = 1.0_FP - step_hermite(x0, x, wx, order) end function real(FP) function get_dphi_shift() !! Finds abscissa, where sigmoid function has value charfun_at_target !! Used to set charfun_at_target correctly real(FP), parameter :: atol = 1.0E-10_FP real(FP) :: x_min, x_max integer :: info type(aux_sigmoid_func_t) :: aux integer, dimension(0) :: iuser ! Dummy array real(FP), dimension(0) :: ruser ! Dummy array if (charfun_at_target < atol) then get_dphi_shift = -0.5_FP * immersed_trace_width_charfun return endif if (charfun_at_target > 1.0_FP - atol) then get_dphi_shift = 0.5_FP * immersed_trace_width_charfun return endif if (immersed_trace_width_charfun < atol) then get_dphi_shift = 0.0_FP return endif x_min = - immersed_trace_width_charfun / 2.0_FP ! Left search boundary x_max = + immersed_trace_width_charfun / 2.0_FP ! Right search boundary call find_zero(x_min, x_max, atol, aux, iuser, ruser, & get_dphi_shift, info) if (info /= 0) then call handle_error('find_zero failed', PARALLAX_ERR_IMMERSED, & __LINE__, __FILE__) endif end function real(FP) function aux_sigmoid_func(this, t, iuser, ruser) !! Auxiliary to get_dphi_shift function class(aux_sigmoid_func_t), intent(in) :: this !! Instance of test function class real(FP), intent(in) :: t ! Abscissa integer, dimension(:), intent(in) :: iuser !! Integer inputs real(FP), dimension(:), intent(in) :: ruser !! Dummy, not used aux_sigmoid_func = 1.0_FP - charfun_at_target & - step_hermite(x0=t, x=0.0_FP, wx=immersed_trace_width_charfun, & order=immersed_step_order) end function subroutine destructor_immersed_trace_t(self) !! Destructor for immersed_rho_t type(immersed_trace_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