submodule (penalisation_m) penalisation_build_s !! Routines related with building penalisation implicit none contains module subroutine build_penalisation(self, comm_handler, equi, mesh, multigrid, dbgout) class(penalisation_t), intent(inout) :: self type(comm_handler_t), intent(in) :: comm_handler class(equilibrium_t) :: equi type(mesh_cart_t), intent(inout) :: mesh type(multigrid_t), intent(inout) :: multigrid integer, intent(in), optional :: dbgout select case(self%pen_method) case(PEN_METHOD_VIA_TRACE) call self%build_via_trace(comm_handler, equi, mesh, dbgout) case(PEN_METHOD_VIA_IN_VESSEL) call self%build_via_in_vessel(comm_handler, equi, mesh, dbgout) case(PEN_METHOD_VIA_RHO) call self%build_via_rho(comm_handler, 'params.in', equi, mesh, dbgout) case(PEN_METHOD_VIA_SUBTRACE) call self%build_via_subtrace(comm_handler, equi, mesh, multigrid, dbgout) case(PEN_METHOD_VIA_STABLE_TRACE) call self%build_via_stable_trace(comm_handler, equi, mesh, dbgout) case(PEN_METHOD_NONE) allocate(self%charfun(mesh%get_n_points_inner()), source=0.0_GP) allocate(self%dirindfun(mesh%get_n_points_inner()), source=0.0_GP) case default call handle_error('Pen_method not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('pen_method', [self%pen_method])) end select ! Build penalisation indices independent of method call self%build_pen_inds(mesh, dbgout) end subroutine module subroutine build_via_trace(self, comm_handler, equi, mesh, dbgout) class(penalisation_t), intent(inout) :: self type(comm_handler_t), intent(in) :: comm_handler class(equilibrium_t) :: equi type(mesh_cart_t), intent(in) :: mesh integer, intent(in), optional :: dbgout real(GP), parameter :: trace_prec = 1.0E-8_GP integer :: outi integer :: ki, l, in_target_i real(GP) :: x, y, rho, phi real(GP) :: dphi_to_fwd_target ! toroidal angle until forward target is hit real(GP) :: dphi_to_bwd_target ! toroidal angle until backward target is hit real(GP) :: dphi_out ! toroidal angle to end of domain (DISTRICT_OUT) real(GP) :: dphi_max ! maximum toroidal angle to be traced real(GP) :: charfun_fwd, charfun_bwd ! characteristic function based on forward/backward target real(GP) :: dphirel real(GP) :: rho_step_char ! for limiter: position of radial step for characteristic function real(GP) :: rho_step_dirind ! for limiter: position of radial step for directional indicator function real(GP) :: dphi_map real(GP) :: dphi_shift ! Shift of input argument for computing characteristic function type(func_1d_charfun_t) :: func_hermite_aux integer, dimension(1) :: iuser real(GP), dimension(2) :: ruser integer :: ifail real(GP) :: x_min, x_max real(GP) :: atol = 1.0E-10_GP ! set output level outi = 0 if (present(dbgout)) then if (is_rank_info_writer) then outi = dbgout endif if (dbgout >= 3) then outi = dbgout ! every rank writes endif endif if (outi > 0) then write(get_stdout(),*)'' write(get_stdout(),*)'Building penalisation VIA_IN_TRACE' endif phi = mesh%get_phi() dphi_map = self%dphi dphi_max = self%dphi_max ! allocate fields allocate(self%charfun(mesh%get_n_points_inner())) allocate(self%dirindfun(mesh%get_n_points_inner())) ! Compute input shift to achieve intended charfun_at_target dphi_shift = GP_NAN select case(self%penfuns_type) case(PENFUNS_TYPE_TANH) ! Set via analytical inverse of tanh stepfunction dphi_shift = - self%charfun_parwidth * atanh(1.0_GP - 2.0_GP * self%charfun_at_target) case(PENFUNS_TYPE_HERMITE) ! No analytic inverse guaranteed, must use numerical search ! Set inputs to wrapper function iuser = [ self%hermite_order ] ruser = [ self%charfun_at_target, self%charfun_parwidth ] x_min = - self%charfun_parwidth / 2.0_GP ! Left search boundary x_max = + self%charfun_parwidth / 2.0_GP ! Right search boundary call find_zero(x_min, x_max, atol, func_hermite_aux, iuser, ruser, dphi_shift, ifail) if ( ifail /= 0 ) then call handle_error('No root within tolerance found', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('x_min, x_max, charfun_at_target, : ', & r_info=[x_min, x_max, self%charfun_at_target])) endif case default call handle_error('penfuns_type not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('penfuns_type: ', [self%penfuns_type])) end select ! Loop over inner grid points do ki = 1, mesh% get_n_points_inner() if (is_rank_info_writer .and. outi >= 1) then call progress_bar(ki, mesh% get_n_points_inner(), 20) endif l = mesh%inner_indices(ki) x = mesh%get_x(l) y = mesh%get_y(l) ! Compute toroidal angles to target ----------------------------------------------------------------------- if (equi%in_vessel(x, y, phi)) then in_target_i = 0 ! Trace until forward / backward target (or maximum distance) ! for inner points phi_to_fwd_target and phi_to_bwd_target are both always positive call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = dphi_max, & equi = equi, & phi_end = dphi_to_fwd_target, & condition = condition_in_target) dphi_to_fwd_target = dphi_to_fwd_target - phi call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = -dphi_max, & equi = equi, & phi_end = dphi_to_bwd_target, & condition = condition_in_target) dphi_to_bwd_target = dphi_to_bwd_target - phi dphi_to_bwd_target = -dphi_to_bwd_target else ! Trace from points within penalisation region in reverse direction until target is hit (or maximum distance) ! for points within penalisation region phi_to_fwd_target and phi_to_bwd_target are both always negative in_target_i = 1 ! For forward target -> trace backwards ------------------------------------------------------------------- ! Trace toroidal distance until field line leaves domain call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = -dphi_max, & equi = equi, & phi_end = dphi_out, & condition = condition_not_in_domain) dphi_out = dphi_out - phi ! Trace until in reverse direction until target is hit (up to at most phi_out) call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = dphi_out, & equi = equi, & phi_end = dphi_to_fwd_target, & condition = condition_not_in_target) dphi_to_fwd_target = dphi_to_fwd_target - phi if ( abs(dphi_out - dphi_to_fwd_target) < trace_prec ) then ! if end of domain is reached before target is hit set distance to target to maximum dphi_to_fwd_target = dphi_max endif ! For backward target -> trace forward, analogously --------------------------------------------------- call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = dphi_max, & equi = equi, & phi_end = dphi_out, & condition = condition_not_in_domain) dphi_out = dphi_out - phi call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = dphi_out, & equi = equi, & phi_end = dphi_to_bwd_target, & condition = condition_not_in_target) dphi_to_bwd_target = dphi_to_bwd_target - phi if ( abs(dphi_out - dphi_to_bwd_target) < trace_prec ) then dphi_to_bwd_target = dphi_max else ! phi_to_bwd_target is defined negative within target dphi_to_bwd_target = -dphi_to_bwd_target endif endif ! Compute characteristic penalisation function ------------------------------------------------------------ select case(self%penfuns_type) case(PENFUNS_TYPE_TANH) ! Charatceristic function based on forward target charfun_fwd = 1.0_GP - step_tanh(x0 = dphi_shift, & x = dphi_to_fwd_target, & wx = self%charfun_parwidth ) ! Charatceristic function based on backward target charfun_bwd = 1.0_GP - step_tanh(x0 = dphi_shift, & x = dphi_to_bwd_target, & wx = self%charfun_parwidth ) case(PENFUNS_TYPE_HERMITE) ! Charatceristic function based on forward target charfun_fwd = 1.0_GP - step_hermite(x0 = dphi_shift, & x = dphi_to_fwd_target, & wx = self%charfun_parwidth, & order = self%hermite_order ) ! Charatceristic function based on backward target charfun_bwd = 1.0_GP - step_hermite(x0 = dphi_shift, & x = dphi_to_bwd_target, & wx = self%charfun_parwidth, & order = self%hermite_order ) case default call handle_error('penfuns_type not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('penfuns_type: ', [self%penfuns_type])) end select self%charfun(ki) = charfun_fwd + charfun_bwd ! Apply correction if forward and backward characteristic functions overlap (e.g. for limiters) if ( (dphi_to_fwd_target < 0.0_GP).and.(dphi_to_bwd_target < 0.0_GP)) then ! Use function that has shortest target to target (maybe instead multiply?) if (dphi_to_fwd_target < dphi_to_bwd_target) then self%charfun(ki) = charfun_bwd else self%charfun(ki) = charfun_fwd endif endif ! Compute direction indication function ------------------------------------------------------------------- dphirel = abs(dphi_to_bwd_target) - abs(dphi_to_fwd_target) select case(self%penfuns_type) case(PENFUNS_TYPE_TANH) self%dirindfun(ki) = -1.0_GP + 2 * step_tanh(x0 = 0.0_GP, & x = dphirel, & wx = self%dirindfun_parwidth) case(PENFUNS_TYPE_HERMITE) self%dirindfun(ki) = -1.0_GP + 2 * step_hermite(x0 = 0.0_GP, & x = dphirel, & wx = self%dirindfun_parwidth, & order = self%hermite_order) case default call handle_error('penfuns_type not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('penfuns_type: ', [self%penfuns_type])) end select enddo ! Radial penalisation for circular limiter geometry select type(equi) type is(circular_t) ! radial step for characteristic function is shifted half width into SOL, i.e. step starts at rho_limiter rho_step_char = equi%rho_limiter + self%charfun_radlimwidth/2.0_GP ! radial step for directional indicator function is at rho_limiter rho_step_dirind = equi%rho_limiter 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%rho(x, y, phi) ! phi = 0 as circular is axisymmetric ! Always hermite step is used due to its compact support self%charfun(ki) = self%charfun(ki) * step_hermite(x0 = rho_step_char, & x = rho, & wx = self%charfun_radlimwidth, & order = self%hermite_order ) ! hard step for direction indicator function self%dirindfun(ki) = self%dirindfun(ki) * heaviside(x0 = rho_step_dirind, x = rho) enddo end select if (outi > 0) then call self%display() endif contains logical function condition_in_target(x, y, phi) real(GP), intent(in) :: x real(GP), intent(in) :: y real(GP), intent(in) :: phi condition_in_target = (.not.equi%in_vessel(x, y, phi)) end function logical function condition_not_in_target(x, y, phi) real(GP), intent(in) :: x real(GP), intent(in) :: y real(GP), intent(in) :: phi condition_not_in_target = equi%in_vessel(x, y, phi) end function logical function condition_not_in_domain(x, y, phi) real(GP), intent(in) :: x real(GP), intent(in) :: y real(GP), intent(in) :: phi integer :: district district = equi%district(x, y, phi) if (district == DISTRICT_OUT) then condition_not_in_domain = .TRUE. else condition_not_in_domain = .FALSE. endif end function end subroutine real(GP) module function func_hermite_aux(this, t, iuser, ruser) class(func_1d_charfun_t), intent(in) :: this real(GP), intent(in) :: t integer, dimension(:), intent(in) :: iuser real(GP), dimension(:), intent(in) :: ruser integer :: hermite_order real(GP) :: charfun_at_target, charfun_parwidth hermite_order = iuser(1) charfun_at_target = ruser(1) charfun_parwidth = ruser(2) func_hermite_aux = 1.0_GP - charfun_at_target & - step_hermite(x0 = t, & x = 0.0_GP, & wx = charfun_parwidth, & order = hermite_order) end function module subroutine build_via_in_vessel(self, comm_handler, equi, mesh, dbgout) class(penalisation_t), intent(inout) :: self type(comm_handler_t), intent(in) :: comm_handler class(equilibrium_t) :: equi type(mesh_cart_t), intent(in) :: mesh integer, intent(in), optional :: dbgout integer :: outi, ki, l real(GP) :: x, y, phi ! set output level outi = 0 if (present(dbgout)) then if (is_rank_info_writer) then outi = dbgout endif if (dbgout >= 3) then outi = dbgout ! every rank writes endif endif if (outi > 0) then write(get_stdout(),*)'' write(get_stdout(),*)'Building penalisation VIA_IN_VESSEL' endif phi = mesh%get_phi() ! allocate fields allocate(self%charfun(mesh%get_n_points_inner())) allocate(self%dirindfun(mesh%get_n_points_inner())) !$omp parallel private(ki, l, x, y) !$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) if (equi%in_vessel(x, y, phi)) then self%charfun(ki) = 0.0_GP ! TODO: Not yet clear what todo with dirindfun self%dirindfun(ki) = 0.0_GP else self%charfun(ki) = 1.0_GP ! TODO: Not yet clear what todo with dirindfun self%dirindfun(ki) = 0.0_GP endif enddo !$omp end do !$omp end parallel end subroutine module subroutine build_via_rho(self, comm_handler, filepath, equi, mesh, dbgout) class(penalisation_t), intent(inout) :: self type(comm_handler_t), intent(in) :: comm_handler character(len=*), intent(in) :: filepath class(equilibrium_t) :: equi type(mesh_cart_t), intent(in) :: mesh integer, intent(in), optional :: dbgout integer :: outi, io_unit, io_error character(len=256) :: io_errmsg integer :: ki, l real(GP) :: phi, rhowidth, x, y, rhon real(GP) :: width_core = -0.1_GP real(GP) :: width_sol = -0.1_GP real(GP) :: step_width_core = 0.0_GP real(GP) :: step_width_sol = 0.0_GP namelist / penalisation_via_rho / & width_core, width_sol, step_width_core, step_width_sol ! Set output level outi = 0 if (present(dbgout)) then if (is_rank_info_writer) then outi = dbgout endif if (dbgout >= 3) then outi = dbgout ! every rank writes endif endif if (outi > 0) then write(get_stdout(),*)'' write(get_stdout(),*)'Building penalisation VIA_RHO' endif ! Read parameters specific for via_rho method from file open(newunit=io_unit, file=filepath, 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 read(io_unit, nml=penalisation_via_rho, iostat=io_error, iomsg=io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, & __LINE__, __FILE__) endif close(io_unit) phi = mesh%get_phi() rhowidth = equi%rhomax - equi%rhomin ! Allocate fields allocate(self%charfun(mesh%get_n_points_inner())) allocate(self%dirindfun(mesh%get_n_points_inner())) !$omp parallel default(none) private(ki, l, x, y, rhon) & !$omp shared(self, equi, mesh, phi, rhowidth, & !$omp width_sol, step_width_sol, width_core, step_width_core) !$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) rhon = (equi%rho(x, y, phi) - equi%rhomin) / rhowidth select case(self%penfuns_type) case(PENFUNS_TYPE_TANH) !use tanh step function self%charfun(ki) = step_tanh(1.0_GP-width_sol, rhon, step_width_sol) & + 1.0_GP - step_tanh(width_core, rhon, step_width_core) case(PENFUNS_TYPE_HERMITE) !use hermite step function self%charfun(ki) = step_hermite(1.0_GP-width_sol, rhon, step_width_sol, order=self%hermite_order) & + 1.0_GP - step_hermite(width_core, rhon, step_width_core, order=self%hermite_order) case default !default to throwing an error message if the function is not recognized call handle_error('penfuns_type not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('penfuns_type: ', [self%penfuns_type])) end select ! Setting dirindfun to 0 is enough for our purpose self%dirindfun(ki) = 0.0_GP enddo !$omp end do !$omp end parallel end subroutine module subroutine build_pen_inds(self, mesh, dbgout) class(penalisation_t), intent(inout) :: self type(mesh_cart_t) :: mesh integer, intent(in), optional :: dbgout integer :: i, ipen, ipenborder, l, ki, kb, kg, outi integer, dimension(mesh%get_n_points_inner()) :: p_inds_buf, & pb_inds_buf integer, dimension(4) :: ind_nb real(GP), dimension(mesh%get_n_points()) :: pen_vals integer, dimension(mesh%get_n_points_boundary()) :: bnd_descrs_nmn ! Set debug output level outi = 0 if ( present(dbgout) ) then if ( is_rank_info_writer ) then outi = dbgout endif if ( dbgout >= 3) then outi = dbgout ! Every rank prints debug info endif endif ipen = 0 ipenborder = 0 ! Extrapolate penalisation values to full mesh !$omp parallel private(l, ki, kb, kg) !$omp do do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) pen_vals(l) = self%get_charfun_val(ki) enddo !$omp end do !$omp do do kb = 1, mesh%get_n_points_boundary() bnd_descrs_nmn(kb) = BND_TYPE_NEUMANN enddo !$omp end do call set_perpbnds(mesh, bnd_descrs_nmn, pen_vals) call extrapolate_ghost_points(mesh, pen_vals) !$omp end parallel do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) if ( pen_vals(l) > GP_EPS ) then ! Find points in penalisation ipen = ipen + 1 p_inds_buf(ipen) = l else ! Find points bordering penalisation call find_neighbor_inds(mesh, l, diagonal=.false., ind_nb=ind_nb) do i = 1, 4 if ( .not. mesh%is_ghost_point(l) ) then ! Exclude ghost points if ( pen_vals(ind_nb(i)) > GP_EPS) then ! Found neighbor point in penalisation region ipenborder = ipenborder + 1 pb_inds_buf(ipenborder) = l endif endif enddo endif enddo ! Allocate and fill actual fields self%n_p_inds = ipen if (.not. allocated(self%p_inds)) allocate(self%p_inds(ipen)) do i = 1, ipen self%p_inds(i) = p_inds_buf(i) enddo self%n_pb_inds = ipenborder if (.not. allocated(self%pb_inds)) allocate(self%pb_inds(ipenborder)) do i = 1, ipenborder self%pb_inds(i) = pb_inds_buf(i) enddo ! Assert that sets are disjunct do i = 1, ipenborder if (any(self%p_inds == self%pb_inds(i))) then call handle_error('Penpoint sets are not disjunct', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('Point in both p_inds & pb_inds: ', & [self%pb_inds(i)])) endif enddo if ( outi > 0 ) then write(get_stdout(),*)'' write(get_stdout(),*)'Info on penalisation indices:' write(get_stdout(),*)'n_p_inds : ', self%n_p_inds write(get_stdout(),*)'n_pb_inds : ', self%n_pb_inds endif end subroutine module subroutine build_via_subtrace(self, comm_handler, equi, mesh, multigrid, dbgout) class(penalisation_t), intent(inout) :: self type(comm_handler_t), intent(in) :: comm_handler class(equilibrium_t) :: equi type(mesh_cart_t), intent(inout) :: mesh type(multigrid_t), intent(inout) :: multigrid integer, intent(in), optional :: dbgout logical :: in_vessel integer :: outi, ntrace, itrace, ki, l, in, jn, ln, sinfo real(GP) :: x, y, xn, yn, res, & dphi_trace, x_trace, y_trace, phi_fwd, phi_trace_fwd, phi_bwd, phi_trace_bwd real(GP), dimension(mesh%get_n_points_inner()) :: wrk_charfun, wrk_dirind, & xfwd, yfwd, xbwd, ybwd real(GP), dimension(mesh%get_n_points_inner()) :: helm_lambda, helm_xi real(GP), dimension(mesh%get_n_points()) :: helm_co, helm_rhs, helm_sol type(parameters_helmholtz_solver_factory) :: helm_params class(helmholtz_solver_t), allocatable :: hsolver ! Set output level outi = 0 if (present(dbgout)) then if (is_rank_info_writer) then outi = dbgout endif if (dbgout >= 3) then outi = dbgout ! every rank writes endif endif if (outi > 0) then write(get_stdout(),*)'' write(get_stdout(),*)'Building penalisation VIA_SUBTRACE' endif ! Build penalisation functions based on subsequent tracing wrk_charfun = 0.0_GP wrk_dirind = 0.0_GP dphi_trace = self%dphi/2.0_GP ntrace = floor(self%dphi_max/dphi_trace) phi_fwd = mesh%get_phi() phi_bwd = mesh%get_phi() do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) xfwd(ki) = mesh%get_x(l) yfwd(ki) = mesh%get_y(l) xbwd(ki) = mesh%get_x(l) ybwd(ki) = mesh%get_y(l) enddo do itrace = 1, ntrace if (outi > 0) then write(get_stdout(),*)' subtracing ', itrace ,' / ', ntrace endif !$omp parallel default(none) & !$omp private(ki, l, x, y, x_trace, y_trace, in_vessel) & !$omp shared(equi, mesh, xfwd, yfwd, xbwd, ybwd, wrk_charfun, wrk_dirind, & !$omp phi_fwd, phi_bwd, dphi_trace, phi_trace_fwd, phi_trace_bwd) !$omp do lastprivate(phi_trace_fwd, phi_trace_bwd) do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) x = mesh%get_x(l) y = mesh%get_y(l) in_vessel = equi%in_vessel(x,y,mesh%get_phi()) call trace(x_init = xfwd(ki), & y_init = yfwd(ki), & phi_init = phi_fwd, & dphi = dphi_trace, & equi = equi, & x_end = x_trace, & y_end = y_trace, & phi_end = phi_trace_fwd) xfwd(ki) = x_trace yfwd(ki) = y_trace phi_trace_fwd = modulo(phi_trace_fwd, TWO_PI) call trace(x_init = xbwd(ki), & y_init = ybwd(ki), & phi_init = phi_bwd, & dphi = -dphi_trace, & equi = equi, & x_end = x_trace, & y_end = y_trace, & phi_end = phi_trace_bwd) xbwd(ki) = x_trace ybwd(ki) = y_trace phi_trace_bwd = modulo(phi_trace_bwd, TWO_PI) ! If point or any traced point goes out of vessel set charfun to 1 if ( (.not.in_vessel) .or. & (.not.equi%in_vessel(xfwd(ki),yfwd(ki),phi_trace_fwd)) .or. & (.not.equi%in_vessel(xbwd(ki),ybwd(ki),phi_trace_bwd)) )then wrk_charfun(ki) = 1.0_GP endif ! If any tracing crosses vessel set dirindfun accordingly if (abs(wrk_dirind(ki)) < GP_EPS) then if (in_vessel .and. (.not.equi%in_vessel(xfwd(ki),yfwd(ki),phi_trace_fwd)) ) then wrk_dirind(ki) = +1.0_GP endif if (in_vessel .and. (.not.equi%in_vessel(xbwd(ki),ybwd(ki),phi_trace_bwd)) ) then wrk_dirind(ki) = -1.0_GP endif if (.not.in_vessel .and. (equi%in_vessel(xfwd(ki),yfwd(ki),phi_trace_fwd)) ) then wrk_dirind(ki) = -1.0_GP endif if (.not.in_vessel .and. (equi%in_vessel(xbwd(ki),ybwd(ki),phi_trace_bwd)) ) then wrk_dirind(ki) = +1.0_GP endif endif enddo !$omp end do !$omp end parallel phi_fwd = phi_trace_fwd phi_bwd = phi_trace_bwd enddo ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !iplane = comm_handler%get_plane() !do ki = 1, mesh%get_n_points_inner() ! write(30+iplane+(cnt-1)*100,*)wrk_charfun(ki),wrk_dirind(ki) !enddo !do ki = 1, mesh%get_n_points_inner() ! read(30+iplane+(cnt-1)*100,*)wrk_charfun(ki),wrk_dirind(ki) !enddo !call MPI_BARRIER(MPI_COMM_WORLD, ierr) ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! Set characteristic function to one for points perpendicular adjacent to vessel do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) if (wrk_charfun(ki) < GP_EPS) then do jn = -mesh%get_size_neighbor(), mesh%get_size_neighbor() do in = -mesh%get_size_neighbor(), mesh%get_size_neighbor() ln = mesh%index_neighbor(in, jn, l) xn = mesh%get_x(ln) yn = mesh%get_y(ln) if (.not.equi%in_vessel(xn,yn,mesh%get_phi())) then wrk_charfun(ki) = 1.0_GP endif enddo enddo endif enddo ! Firstly build penalisation functions via in_vessel call self%build_via_in_vessel(comm_handler, equi, mesh, dbgout=0) ! Interpolate region between vessel and region where wrk_charfun>0 via helmholtz solve ! ... for charfun helm_co = 1.0_GP helm_rhs = 0.0_GP helm_sol = 0.0_GP do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) if ( (self%charfun(ki) > 0.9_GP) .or. (wrk_charfun(ki) < 0.1_GP) ) then helm_lambda(ki) = 1.0_GP helm_xi(ki) = 0.0_GP helm_rhs(l) = self%charfun(ki) else helm_lambda(ki) = 0.0_GP helm_xi(ki) = 1.0_GP helm_rhs(l) = 0.0_GP endif enddo call helmholtz_solver_factory(multigrid, & BND_TYPE_NEUMANN, & BND_TYPE_NEUMANN, & BND_TYPE_NEUMANN, & BND_TYPE_NEUMANN, & helm_co, helm_lambda, helm_xi, & helm_params, 'DIRECT', hsolver) call hsolver%solve(helm_rhs, helm_sol, res, sinfo) deallocate(hsolver) if (sinfo < 0) then call handle_error('Solver for penalisation failed', & GRILLIX_ERR_SOLVER2D, __LINE__, __FILE__, & error_info_t('sinfo, res: ',[sinfo],[res])) endif !$omp parallel default(none) private(ki,l) shared(self, mesh, helm_sol) !$omp do do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) self%charfun(ki) = helm_sol(l) enddo !$omp end do !$omp end parallel ! Interpolate dirindfun via helmholtz solve helm_co = 1.0_GP helm_rhs = 0.0_GP helm_sol = 0.0_GP do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) if (abs(wrk_dirind(ki)) > 0.1_GP) then helm_lambda(ki) = 1.0_GP helm_xi(ki) = 0.0_GP helm_rhs(l) = wrk_dirind(ki) else helm_lambda(ki) = 0.0_GP helm_xi(ki) = 1.0_GP helm_rhs(l) = 0.0_GP endif enddo call helmholtz_solver_factory(multigrid, & BND_TYPE_NEUMANN, & BND_TYPE_NEUMANN, & BND_TYPE_NEUMANN, & BND_TYPE_NEUMANN, & helm_co, helm_lambda, helm_xi, & helm_params, 'DIRECT', hsolver) call hsolver%solve(helm_rhs, helm_sol, res, sinfo) deallocate(hsolver) if (sinfo < 0) then call handle_error('Solver for penalisation failed', & GRILLIX_ERR_SOLVER2D, __LINE__, __FILE__, & error_info_t('sinfo, res: ',[sinfo],[res])) endif !$omp parallel default(none) private(ki,l) shared(self, mesh, helm_sol) !$omp do do ki = 1, mesh%get_n_points_inner() l = mesh%inner_indices(ki) self%dirindfun(ki) = helm_sol(l) enddo !$omp end do !$omp end parallel if (outi > 0) then call self%display() endif end subroutine module subroutine destructor(self) type(penalisation_t), intent(inout) :: self self%pen_method = 0 self%penfuns_type = 0 self%hermite_order = 0 self%charfun_parwidth = GP_NAN self%dirindfun_parwidth = GP_NAN if (allocated(self%dirindfun)) deallocate(self%dirindfun) if (allocated(self%charfun)) deallocate(self%charfun) if (allocated(self%p_inds)) deallocate(self%p_inds) if (allocated(self%pb_inds)) deallocate(self%pb_inds) end subroutine module subroutine build_via_stable_trace(self, comm_handler, equi, mesh, dbgout) class(penalisation_t), intent(inout) :: self type(comm_handler_t), intent(in) :: comm_handler class(equilibrium_t) :: equi type(mesh_cart_t), intent(in) :: mesh integer, intent(in), optional :: dbgout real(GP), parameter :: trace_prec = 1.0E-8_GP integer :: outi integer :: ki, l, in_target_i real(GP) :: x, y, rho, phi real(GP) :: dphi_to_fwd_target ! toroidal angle until forward target is hit real(GP) :: dphi_to_bwd_target ! toroidal angle until backward target is hit real(GP) :: dphi_out ! toroidal angle to end of domain (DISTRICT_OUT) real(GP) :: dphi_max ! maximum toroidal angle to be traced real(GP) :: charfun_fwd, charfun_bwd ! characteristic function based on forward/backward target real(GP) :: dphirel real(GP) :: rho_step_char ! for limiter: position of radial step for characteristic function real(GP) :: rho_step_dirind ! for limiter: position of radial step for directional indicator function real(GP) :: dphi_map real(GP) :: dphi_shift ! Shift of input argument for computing characteristic function type(func_1d_charfun_t) :: func_hermite_aux integer, dimension(1) :: iuser real(GP), dimension(2) :: ruser integer :: ifail real(GP) :: x_min, x_max real(GP) :: atol = 1.0E-10_GP ! set output level outi = 0 if (present(dbgout)) then if (is_rank_info_writer) then outi = dbgout endif if (dbgout >= 3) then outi = dbgout ! every rank writes endif endif if (outi > 0) then write(get_stdout(),*)'' write(get_stdout(),*)'Building penalisation VIA_STABLE_TRACE' endif phi = mesh%get_phi() dphi_map = self%dphi dphi_max = self%dphi_max ! allocate fields allocate(self%charfun(mesh%get_n_points_inner())) allocate(self%dirindfun(mesh%get_n_points_inner())) ! Compute input shift to achieve intended charfun_at_target dphi_shift = GP_NAN select case(self%penfuns_type) case(PENFUNS_TYPE_HERMITE) ! No analytic inverse guaranteed, must use numerical search ! Set inputs to wrapper function iuser = [ self%hermite_order ] ruser = [ self%charfun_at_target, self%charfun_parwidth ] x_min = - self%charfun_parwidth / 2.0_GP ! Left search boundary x_max = + self%charfun_parwidth / 2.0_GP ! Right search boundary call find_zero(x_min, x_max, atol, func_hermite_aux, iuser, ruser, dphi_shift, ifail) if ( ifail /= 0 ) then call handle_error('No root within tolerance found', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('x_min, x_max, charfun_at_target, : ', & r_info=[x_min, x_max, self%charfun_at_target])) endif case default call handle_error('penfuns_type not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('penfuns_type: ', [self%penfuns_type])) end select ! Loop over inner grid points do ki = 1, mesh% get_n_points_inner() if (is_rank_info_writer .and. outi >= 1) then call progress_bar(ki, mesh% get_n_points_inner(), 20) endif l = mesh%inner_indices(ki) x = mesh%get_x(l) y = mesh%get_y(l) ! Compute toroidal angles to target ----------------------------------------------------------------------- if (equi%in_vessel(x, y, phi)) then if(equi%rho(x, y, phi) < self%rho_min) then self%charfun(ki) = 0.0_GP self%dirindfun(ki) = 0.0_GP cycle endif in_target_i = 0 ! Trace until forward / backward target (or maximum distance) ! for inner points phi_to_fwd_target and phi_to_bwd_target are both always positive call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = dphi_max, & equi = equi, & phi_end = dphi_to_fwd_target, & condition = condition_in_target, & stop_int = .true., & maxstep = self%max_step_size) dphi_to_fwd_target = dphi_to_fwd_target - phi call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = -dphi_max, & equi = equi, & phi_end = dphi_to_bwd_target, & condition = condition_in_target, & stop_int = .true., & maxstep = self%max_step_size) dphi_to_bwd_target = dphi_to_bwd_target - phi dphi_to_bwd_target = -dphi_to_bwd_target else ! Trace from points within penalisation region in reverse direction until target is hit (or maximum distance) ! for points within penalisation region phi_to_fwd_target and phi_to_bwd_target are both always negative in_target_i = 1 ! For forward target -> trace backwards ------------------------------------------------------------------- ! Trace toroidal distance until field line leaves domain call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = -dphi_max, & equi = equi, & phi_end = dphi_out, & condition = condition_not_in_domain, & stop_int = .true., & maxstep = self%max_step_size) dphi_out = dphi_out - phi ! Trace until in reverse direction until target is hit (up to at most phi_out) call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = dphi_out, & equi = equi, & phi_end = dphi_to_fwd_target, & condition = condition_not_in_target, & stop_int = .true., & maxstep = self%max_step_size) dphi_to_fwd_target = dphi_to_fwd_target - phi if ( abs(dphi_out - dphi_to_fwd_target) < trace_prec ) then ! if end of domain is reached before target is hit set distance to target to maximum dphi_to_fwd_target = dphi_max endif ! For backward target -> trace forward, analogously --------------------------------------------------- call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = dphi_max, & equi = equi, & phi_end = dphi_out, & condition = condition_not_in_domain, & stop_int = .true., & maxstep = self%max_step_size) dphi_out = dphi_out - phi call trace(x_init = x, & y_init = y, & phi_init = phi, & dphi = dphi_out, & equi = equi, & phi_end = dphi_to_bwd_target, & condition = condition_not_in_target, & stop_int = .true., & maxstep = self%max_step_size) dphi_to_bwd_target = dphi_to_bwd_target - phi if ( abs(dphi_out - dphi_to_bwd_target) < trace_prec ) then dphi_to_bwd_target = dphi_max else ! phi_to_bwd_target is defined negative within target dphi_to_bwd_target = -dphi_to_bwd_target endif endif ! Compute characteristic penalisation function ------------------------------------------------------------ select case(self%penfuns_type) case(PENFUNS_TYPE_HERMITE) ! Charatceristic function based on forward target charfun_fwd = 1.0_GP - step_hermite(x0 = dphi_shift, & x = dphi_to_fwd_target, & wx = self%charfun_parwidth, & order = self%hermite_order ) ! Charatceristic function based on backward target charfun_bwd = 1.0_GP - step_hermite(x0 = dphi_shift, & x = dphi_to_bwd_target, & wx = self%charfun_parwidth, & order = self%hermite_order ) case default call handle_error('penfuns_type not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('penfuns_type: ', [self%penfuns_type])) end select self%charfun(ki) = charfun_fwd + charfun_bwd ! Apply correction if forward and backward characteristic functions overlap (e.g. for limiters) if ( (dphi_to_fwd_target < 0.0_GP).and.(dphi_to_bwd_target < 0.0_GP)) then ! Use function that has shortest target to target (maybe instead multiply?) if (dphi_to_fwd_target < dphi_to_bwd_target) then self%charfun(ki) = charfun_bwd else self%charfun(ki) = charfun_fwd endif endif ! Compute direction indication function ------------------------------------------------------------------- dphirel = abs(dphi_to_bwd_target) - abs(dphi_to_fwd_target) select case(self%penfuns_type) case(PENFUNS_TYPE_HERMITE) self%dirindfun(ki) = -1.0_GP + 2 * step_hermite(x0 = 0.0_GP, & x = dphirel, & wx = self%dirindfun_parwidth, & order = self%hermite_order) case default call handle_error('penfuns_type not valid', & GRILLIX_ERR_OTHER, __LINE__, __FILE__, & error_info_t('penfuns_type: ', [self%penfuns_type])) end select enddo if (outi > 0) then call self%display() endif contains logical function condition_in_target(x, y, phi) real(GP), intent(in) :: x real(GP), intent(in) :: y real(GP), intent(in) :: phi condition_in_target = (.not.equi%in_vessel(x, y, phi)) end function logical function condition_not_in_target(x, y, phi) real(GP), intent(in) :: x real(GP), intent(in) :: y real(GP), intent(in) :: phi condition_not_in_target = equi%in_vessel(x, y, phi) end function logical function condition_not_in_domain(x, y, phi) real(GP), intent(in) :: x real(GP), intent(in) :: y real(GP), intent(in) :: phi integer :: district district = equi%district(x, y, phi) if (district == DISTRICT_OUT) then condition_not_in_domain = .TRUE. else condition_not_in_domain = .FALSE. endif end function end subroutine end submodule