module parbnd_taylor_m !! Handles setting of parallel boundary conditions !! according to Taylor expansion along boundary conditions !! See Appendix B of !! Andreas Stegmeir et al. 2018 Plasma Phys. Control. Fusion 60 035005. !! DOI 10.1088/1361-6587/aaa373 use precision_m, only : FP, FP_LARGEST, FP_NAN, FP_EPS use screen_io_m, only : get_stdout use constants_m, only : PI use equilibrium_m, only : equilibrium_t use mesh_cart_m, only : mesh_cart_t use connection_length_m, only : distance_to_boundary use fieldline_tracer_m, only : trace use error_handling_m, only : handle_error use status_codes_m, only : PARALLAX_ERR_PARBND, & PARALLAX_ERR_PARAMETERS use elementary_functions_m, only : step_hermite implicit none private integer, parameter, private :: max_taylor_order = 3 real(FP), private :: dphi_max = PI / 2.0_FP !! Sets depth of guard points layer real(FP), protected :: trace_max_step_size = PI / 180.0_FP !! Maximum step size of trace algorithm real(FP), protected :: depth_alt = 0.5_FP !! Depth at which method switches to alternate expression real(FP), protected :: wd_alt = 0.5_FP !! Width over which method switches to alternate expression integer, protected :: hermite_order_alt = 2 !! Order of transition (hermite polynomials) to alternate expression real(FP), protected :: depth_flb = 1.5_FP !! Depth at which method switches to fallback expression real(FP), protected :: wd_flb = 0.5_FP !! Width over which method switches to fallback expression integer, protected :: hermite_order_flb = 2 !! Order of transition (hermite polynomials) to fallback expression real(FP), protected :: depth_ctgy = 2.5_FP !! Depth at which method switches to contingency (deep depth) expression real(FP), protected :: wd_ctgy = 0.5_FP !! Width over which method switches to contingency expression integer, protected :: hermite_order_ctgy = 2 !! Order of transition (hermite polynomials) to contingency expression namelist / parbnd_taylor / & dphi_max, & trace_max_step_size, & depth_alt, wd_alt, hermite_order_alt, & depth_flb, wd_flb, hermite_order_flb, & depth_ctgy, wd_ctgy, hermite_order_ctgy type, public :: parbnd_taylor_t !! Handles setting of boundary conditions at parallel boundaries !! according to Taylor expansion method integer, public, allocatable, dimension(:) :: b_bnd_direction !! Indicates if the magnetic field is directed towards or away the !! boundary at the closest parallel boundary: !! - = 1 : if B points toward the boundary !! - = -1 : if B points away from the boundary) !! - = 0 : if no boundary is found nearby (e.g., in closed field line !! region or when the boundary is too far from the point) real(FP), public, allocatable, dimension(:,:) :: dists_to_bnd !! Distances along the magnetic field line to the boundary !! !! For each point l: !! - dists_to_bnd(l, 1): distance from point l to the boundary !! - dists_to_bnd(l, 2): distance from the mapped point of point l !! to the boundary !! - dists_to_bnd(l, 3): distance from the next-but-one mapped !! point of point l to the boundary !! !! Sign convention: !! - Distance is positive if measured in the direction of the magnetic !! field B !! - Distance is negative if measured in the opposite direction of B real(FP), public, allocatable, dimension(:) :: depth_in_bnd !! Corresponding toridal angle to boundary in units of dphi !! - Positive values indicate that point is inside boundary !! - Negative values indicate that point is outside boundary !! - A value of zero indicates that boundary is too far away !! (e.g., in closed field line region) contains procedure, public :: init => init_parbnd_taylor procedure, private :: read_params_parbnd_taylor procedure, public :: display procedure, public :: write_netcdf => write_netcdf_parbnd_taylor procedure, public :: read_netcdf => read_netcdf_parbnd_taylor procedure, private :: taylor_determinant procedure, private :: alt_taylor_determinant procedure, public :: set_guard_points procedure, public :: extrapolate_guard_points final :: destructor_parbnd_taylor_t end type interface module subroutine write_netcdf_parbnd_taylor(self, fgid) !! Writes parbnd_taylor data to NetCDF id class(parbnd_taylor_t), intent(in) :: self !! Instance of class integer, intent(in) :: fgid !! NetCDF file or group id end subroutine module subroutine read_netcdf_parbnd_taylor(self, fgid) !! Reads parbnd_taylor data from NetCDF id class(parbnd_taylor_t), intent(inout) :: self !! Instance of class integer, intent(in) :: fgid !! NetCDF file or group id end subroutine end interface contains subroutine init_parbnd_taylor(self, equi, mesh, dphi, filename) !! Initialises parbnd_taylor, by computing distances to boundaries class(parbnd_taylor_t), intent(inout) :: self !! Instance of type class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh !! Mesh real(FP), intent(in) :: dphi !! Toroidal grid distance character(len=*), intent(in), optional :: filename !! Filename, where to read parameters from, !! if not provided default parameters will be used integer :: l real(FP) :: x, y, phi, arclength_fwd, arclength_bwd, dphi_fwd, & dphi_bwd, arclen_map, dphi_map, x_upstream, y_upstream, phi_upstream logical :: target_reached_fwd, target_reached_bwd if (present(filename)) then call self%read_params_parbnd_taylor(filename) endif allocate(self%b_bnd_direction(mesh%get_n_points())) allocate(self%dists_to_bnd(mesh%get_n_points(), max_taylor_order)) allocate(self%depth_in_bnd(mesh%get_n_points())) phi = mesh%get_phi() !$omp parallel default(none) & !$omp private(l, x, y, arclength_fwd, arclength_bwd, & !$omp target_reached_fwd, target_reached_bwd, dphi_fwd, & !$omp dphi_bwd, dphi_map, x_upstream, y_upstream, phi_upstream, & !$omp arclen_map) & !$omp shared(self, equi, mesh, phi, dphi, dphi_max, trace_max_step_size) !$omp do do l = 1, mesh%get_n_points() x = mesh%get_x(l) y = mesh%get_y(l) self%b_bnd_direction(l) = 0 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) if ((.not.target_reached_fwd) .and. (.not.target_reached_bwd)) then self%dists_to_bnd(l,:) = sqrt(FP_LARGEST) self%depth_in_bnd(l) = 0.0_FP cycle endif ! Determine if B points toward or away from closest boundary if ((target_reached_fwd) .and. (.not.target_reached_bwd)) then self%b_bnd_direction(l) = 1 elseif ((target_reached_bwd) .and. (.not.target_reached_fwd)) then self%b_bnd_direction(l) = -1 else if (abs(arclength_fwd) <= abs(arclength_bwd)) then self%b_bnd_direction(l) = 1 else self%b_bnd_direction(l) = -1 endif endif ! Determine depth in boundary if (self%b_bnd_direction(l) == 1) then self%depth_in_bnd(l) = -dphi_fwd / dphi elseif (self%b_bnd_direction(l) == -1) then self%depth_in_bnd(l) = -dphi_bwd / dphi endif ! Distance to boundary with above described sign convention if (self%b_bnd_direction(l) == 1) then self%dists_to_bnd(l, 1) = -arclength_fwd elseif (self%b_bnd_direction(l) == -1) then self%dists_to_bnd(l, 1) = arclength_bwd endif ! Trace to upstream map point ! and compute its distance to boundary according to sign convention dphi_map = -self%b_bnd_direction(l) * dphi call trace(x, y, phi, dphi_map, equi, & x_end=x_upstream, y_end=y_upstream, phi_end=phi_upstream, & arclen=arclen_map) self%dists_to_bnd(l, 2) = self%dists_to_bnd(l, 1) & - self%b_bnd_direction(l) * arclen_map ! Trace further to next-but-one upstream map point ! and compute its distance to boundary according to sign convention call trace(x_upstream, y_upstream, phi_upstream, dphi_map, equi, & arclen=arclen_map) self%dists_to_bnd(l, 3) = self%dists_to_bnd(l, 2) & - self%b_bnd_direction(l) * arclen_map enddo !$omp end do !$omp end parallel end subroutine subroutine read_params_parbnd_taylor(self, filename) !! Reads parameters for parbnd_taylor class(parbnd_taylor_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=parbnd_taylor, 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) !! Writes basic information of parbnd_taylor_t to stdout class(parbnd_taylor_t), intent(in) :: self !! Instance of class integer :: io_unit, io_error character(len=256) :: io_errmsg io_unit = get_stdout() write(io_unit, nml=parbnd_taylor, 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 set_guard_points(self, equi, mesh, fac_dir, fac_nmn, bnd_vals, & u, u_bwd2, u_bwd1, u_fwd1, u_fwd2, order, uisect) !! Sets boundary condition on guard points according to eq. (B.2) !! in DOI 10.1088/1361-6587/aaa373, with modifications for deep guard !! points class(parbnd_taylor_t), intent(in) :: self !! Instance of class class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh !! Mesh real(FP), intent(in) :: fac_dir !! Factor for Dirichlet boundary condition !! (beta in DOI 10.1088/1361-6587/aaa373) real(FP), intent(in) :: fac_nmn !! Factor for Neumann boundary condition !! (gamma in DOI 10.1088/1361-6587/aaa373) real(FP), intent(in), dimension(mesh%get_n_points()) :: bnd_vals !! Boundary values at intersection point of field line with boundary !! (alpha in DOI 10.1088/1361-6587/aaa373) real(FP), dimension(mesh%get_n_points()), intent(inout) :: u !! Field on which guard points are set according to boundary condition real(FP), dimension(mesh%get_n_points()), intent(in) :: & u_bwd2, u_bwd1, u_fwd1, u_fwd2 !! Map points up to two backward and two forward planes integer, intent(in), optional :: order !! Specifies the order of the method (1 or 2, default 2) real(FP), dimension(mesh%get_n_points()), intent(out), optional :: & uisect !! Values of field at intersection point of field line with boundary integer :: order_set, l real(FP) :: x, y, phi, ua, ub, uc, ux, sa, sb, sc, det, det_alt, & ua_alt, ux_alt, ua_flb, ua_ctgy, env_alt, env_flb, env_ctgy, & ua_star, ua_dag, bval logical :: in_vessel if (present(order)) then order_set = order else order_set = 2 endif phi = mesh%get_phi() !$omp parallel default(none) & !$omp private(l, x, y, in_vessel, bval, ua, ub, uc, ux, sa, sb, sc, & !$omp det, det_alt, ua_alt, ux_alt, ua_flb, ua_ctgy, & !$omp env_alt, env_flb, env_ctgy, ua_star, ua_dag) & !$omp shared(self, equi, mesh, phi, bnd_vals, fac_dir, fac_nmn, & !$omp order_set, depth_alt, wd_alt, hermite_order_alt, & !$omp depth_flb, wd_flb, hermite_order_flb, & !$omp depth_ctgy, wd_ctgy, hermite_order_ctgy, & !$omp u_bwd1, u_bwd2, u_fwd1, u_fwd2, u, uisect) !$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) if (in_vessel) then if (present(uisect)) then uisect(l) = 0.0_FP endif cycle endif if (self%b_bnd_direction(l) == 0) then if (present(uisect)) then uisect(l) = 0.0_FP endif cycle elseif (self%b_bnd_direction(l) == 1) then ub = u_bwd1(l) uc = u_bwd2(l) elseif (self%b_bnd_direction(l) == -1) then ub = u_fwd1(l) uc = u_fwd2(l) else !$omp critical call handle_error('Value for b_bnd_direction not valid', & PARALLAX_ERR_PARBND, & __LINE__, __FILE__) !$omp end critical endif bval = bnd_vals(l) sa = self%dists_to_bnd(l,1) sb = self%dists_to_bnd(l,2) sc = self%dists_to_bnd(l,3) det = self%taylor_determinant(l, fac_dir, fac_nmn, order_set) if (order_set == 1) then ua = (sb - sa) * (bval) + (fac_dir * sa - fac_nmn) * ub ua = ua / det ux = (bval * sb - fac_nmn * ub) / det elseif (order_set == 2) then ua = (sa - sb) * (sa - sc) * (sc - sb) * bval & - (fac_dir * sa * sc - fac_nmn * (sa + sc)) & * (sa - sc) * ub & + (fac_dir * sa * sb - fac_nmn * (sa + sb)) & * (sa - sb) * uc ua = ua / (2.0_FP * det) ux = sb * sc * (sc - sb) * bval & - fac_nmn * (sc**2 * ub - sb**2 * uc) ux = ux / (2.0_FP * det) endif ! Alternate value, first order skipping adjacent plane det_alt = self%alt_taylor_determinant(l, fac_dir, fac_nmn) ua_alt = (sc - sa) * (bval) + (fac_dir * sa - fac_nmn) * uc ua_alt = ua_alt / det_alt ux_alt = (bval * sc - fac_nmn * uc) / det_alt ! Linear extrapolation as fallback for deeper guard point region ua_flb = 2.0_FP * ub - uc ! Use upstream value as contingency for even deeper guard region ua_ctgy = ub ! Envelope functions env_alt = step_hermite(depth_alt, abs(self%depth_in_bnd(l)), & wd_alt, hermite_order_alt) env_flb = step_hermite(depth_flb, abs(self%depth_in_bnd(l)), & wd_flb, hermite_order_flb) env_ctgy = step_hermite(depth_ctgy, abs(self%depth_in_bnd(l)), & wd_ctgy, hermite_order_ctgy) ! Final expression ua_star = (1.0_FP - env_alt) * ua + env_alt * ua_alt ua_dag = (1.0_FP - env_flb) * ua_star + env_flb * ua_flb u(l) = (1.0_FP - env_ctgy) * ua_dag + env_ctgy * ua_ctgy if (present(uisect)) then uisect(l) = (1.0_FP - env_alt) * ux + env_alt * ux_alt endif enddo !$omp end do !$omp end parallel end subroutine subroutine extrapolate_guard_points(self, equi, mesh, & u, u_bwd2, u_bwd1, u_fwd1, u_fwd2, uisect) !! Extrapolates field into parallel guard points region. !! Method is based on prescribing second parallel derivative along field !! lines to zero. class(parbnd_taylor_t), intent(in) :: self !! Instance of class class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh !! Mesh real(FP), dimension(mesh%get_n_points()), intent(inout) :: u !! Field on which guard points are set according to boundary condition real(FP), dimension(mesh%get_n_points()), intent(in) :: & u_bwd2, u_bwd1, u_fwd1, u_fwd2 !! Map points up to two backward and two forward planes real(FP), dimension(mesh%get_n_points()), intent(out), optional :: & uisect !! Values of field at intersection point of field line with boundary logical :: in_vessel integer :: l real(FP) :: x, y, phi, ua, ub, uc, ux, sa, sb, sc, den, ua_ctgy, & env_ctgy phi = mesh%get_phi() !$omp parallel default(none) & !$omp private(l, x, y, in_vessel, ua, ub, uc, ux, sa, sb, sc, den, & !$omp ua_ctgy, env_ctgy) & !$omp shared(self, equi, mesh, phi, u_bwd1, u_bwd2, u_fwd1, u_fwd2, u, & !$omp uisect, depth_ctgy, wd_ctgy, hermite_order_ctgy) !$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) if (in_vessel) then if (present(uisect)) then uisect(l) = 0.0_FP endif cycle endif if (self%b_bnd_direction(l) == 0) then if (present(uisect)) then uisect(l) = 0.0_FP endif cycle elseif (self%b_bnd_direction(l) == 1) then ub = u_bwd1(l) uc = u_bwd2(l) elseif (self%b_bnd_direction(l) == -1) then ub = u_fwd1(l) uc = u_fwd2(l) else !$omp critical call handle_error('Value for b_bnd_direction not valid', & PARALLAX_ERR_PARBND, & __LINE__, __FILE__) !$omp end critical endif sa = self%dists_to_bnd(l,1) sb = self%dists_to_bnd(l,2) sc = self%dists_to_bnd(l,3) ! Value based on Taylor expansion with zero second derivative den = 1.0_FP / (sb - sc) ua = ((sa - sc) * ub - (sa - sb) * uc) * den ! Use upstream value as contingency for very deep guard region ua_ctgy = ub ! Envelope functions env_ctgy = step_hermite(depth_ctgy, abs(self%depth_in_bnd(l)), & wd_ctgy, hermite_order_ctgy) ! Final expression u(l) = (1.0_FP - env_ctgy) * ua + env_ctgy * ua_ctgy ux = (-sc * ub + sb * uc) * den if (present(uisect)) then uisect(l) = ux endif enddo !$omp end do !$omp end parallel end subroutine pure real(FP) function taylor_determinant(self, l, fac_dir, fac_nmn, order) !! Computes determinant of underlying linear equation system class(parbnd_taylor_t), intent(in) :: self !! Instance of class integer, intent(in) :: l !! Mesh index real(FP), intent(in) :: fac_dir !! Factor for Dirichlet boundary condition real(FP), intent(in) :: fac_nmn !! Factor for Neumann boundary condition integer, intent(in) :: order !! Specifies the order of the method (allowed values: 1 or 2). !! If a value other than 1 or 2 is provided, the function returns NaN. real(FP) :: sb, sc sb = self%dists_to_bnd(l, 2) sc = self%dists_to_bnd(l, 3) if (order == 1) then taylor_determinant = fac_dir * sb - fac_nmn elseif (order == 2) then taylor_determinant = (fac_dir * sb * sc - fac_nmn * (sb + sc)) taylor_determinant = 0.5_FP * taylor_determinant * (sc - sb) else taylor_determinant = FP_NAN endif end function pure real(FP) function alt_taylor_determinant(self, l, fac_dir, fac_nmn) !! Computes determinant of underlying linear equation system !! for alternate expression, i.e. first order Taylor expansion skipping !! adjacent plane class(parbnd_taylor_t), intent(in) :: self !! Instance of class integer, intent(in) :: l !! Mesh index real(FP), intent(in) :: fac_dir !! Factor for Dirichlet boundary condition real(FP), intent(in) :: fac_nmn !! Factor for Neumann boundary condition real(FP) :: sc sc = self%dists_to_bnd(l, 3) alt_taylor_determinant = fac_dir * sc - fac_nmn end function subroutine destructor_parbnd_taylor_t(self) !! Destructor for parbnd_taylor_t type(parbnd_taylor_t), intent(inout) :: self if (allocated(self%b_bnd_direction)) deallocate(self%b_bnd_direction) if (allocated(self%dists_to_bnd)) deallocate(self%dists_to_bnd) if (allocated(self%depth_in_bnd)) deallocate(self%depth_in_bnd) end subroutine end module