module fieldline_tracer_m !! Field line tracing via the DOP853 integrator use screen_io_m, only: get_stderr use comm_handling_m, only: is_master use precision_m, only: FP, FP_EPS, DP, DP_EPS use error_handling_m, only: handle_error, error_info_t use status_codes_m, only: PARALLAX_SUCCESS, PARALLAX_WRN_TRACE, & PARALLAX_ERR_TRACE, PARALLAX_ERR_PARAMETERS use equilibrium_m, only: equilibrium_t use dop853_module use dop853_constants implicit none private public :: trace public :: set_tolerance integer, parameter :: neq = 5 real(DP), save :: rtol(neq) = & [1.0E-8_DP, DP_EPS, DP_EPS, 1.0E-8_DP, 1.0E-8_DP] !! Relative error tolerance for DOP853. Tracing for x and y coordinates !! restricted by absolute error, lengths by relative error. real(DP), save :: atol(neq) = & [DP_EPS, 1.0E-8_DP, 1.0E-8_DP, DP_EPS, DP_EPS] !! Absolute error tolerance for DOP853. Tracing for x and y coordinates !! restricted by absolute error, lengths by relative error. type, extends(dop853_class) :: dop853_equi_t !! DOP integrator which can trace the magnetic field lines of a given !! equilibrium class(equilibrium_t), private, pointer :: equi => null() procedure(condition_interface), private, pointer, nopass :: & condition => condition_default logical, public :: stop_integration = .false. !! Control method how to stop integration, when condition is met. !! If = true, then the routine stops the integration once the !! condition is met. Otherwise, integration continues with integrand !! multiplied by zero. The default setting is false. contains procedure :: set_equilibrium procedure :: set_condition procedure :: fcn => integrands procedure :: solout => stop_tracing end type dop853_equi_t abstract interface logical function condition_interface(x, y, phi) import FP real(FP), intent(in) :: x, y, phi end function condition_interface end interface contains logical function condition_default(x, y, phi) real(FP), intent(in) :: x, y, phi condition_default =.false. end function condition_default subroutine set_equilibrium(this, equi) class(dop853_equi_t), intent(inout) :: this class(equilibrium_t), target, intent(in) :: equi this%equi => equi end subroutine set_equilibrium subroutine set_condition(this, condition) class(dop853_equi_t), intent(inout) :: this procedure(condition_interface) :: condition this%condition => condition end subroutine set_condition subroutine integrands(me, x, y, f) !! Integrand of field line tracing for DOP853 routine class(dop853_equi_t), intent(inout) :: me !! Instance of the type real(DP), intent(in) :: x !! Toroidal angle phi real(DP), dimension(:), intent(in) :: y !! y(1): current phi !! y(2): current x-coordinate !! y(3): current y-coordinate !! y(4): current arclength along field line !! y(5): current 1 / Btor integrated along field line real(DP), dimension(:), intent(out) :: f !! f(1): 1 !! f(2): bx / btor !! f(3): by / btor !! f(4): abs(B / Btor) * Jacobian, !! f(5): 1 / Btor * Jacobian !! Set to zero if the trace condition is true real(DP) :: jac, bx, by, btor, fac real(FP) :: xfp, yfp, phifp phifp = real(x, kind=FP) xfp = real(y(2), kind=FP) yfp = real(y(3), kind=FP) if (me%condition(xfp, yfp, phifp) .and. .not. me%stop_integration) then fac = 0.0_DP else fac = 1.0_DP endif bx = real(me%equi%bx(xfp, yfp, phifp), kind=DP) by = real(me%equi%by(xfp, yfp, phifp), kind=DP) btor = real(me%equi%btor(xfp, yfp, phifp), kind=DP) jac = real(me%equi%jacobian(xfp, yfp, phifp), kind=DP) f(1) = fac f(2) = bx / btor * jac * fac f(3) = by / btor * jac * fac f(4) = sqrt( bx * bx / (btor * btor) + by * by / (btor * btor) & + 1.0_DP) * jac * fac f(5) = 1.0_DP / btor * jac * fac end subroutine integrands subroutine stop_tracing(me, nr, xold, x, y, irtrn, xout) !! A routine to stop the integration once the condition is fulfilled class(dop853_equi_t),intent(inout) :: me !! Instance of the type integer,intent(in) :: nr !! grid point (0,1,...) real(DP),intent(in) :: xold !! the preceding grid point real(DP),intent(in) :: x !! Toroidal angle phi real(DP),dimension(:),intent(in) :: y !! y(1): current phi !! y(2): current x-coordinate !! y(3): current y-coordinate !! y(4): current arclength along field line !! y(5): current 1 / Btor integrated along field line integer,intent(inout) :: irtrn !! Serves to interrupt the integration. if `irtrn` is set `<0`, !! [[dop853]] will return to the calling program. real(DP),intent(out) :: xout !! Can be used for efficient intermediate output if (me%condition(y(2), y(3), x)) then irtrn = -1 endif end subroutine stop_tracing subroutine set_tolerance(rtol_, atol_) !! Sets the DOP853 relative and absolute error tolerances for each !! equation solved during the field line tracing real(DP), dimension(neq), intent(in) :: rtol_ !! Relative error tolerance for each equation real(DP), dimension(neq), intent(in) :: atol_ !! Absolute error tolerance for each equation rtol = rtol_ atol = atol_ end subroutine set_tolerance subroutine trace(x_init, y_init, phi_init, dphi, equi, x_end, y_end, & phi_end, arclen, fluxexp, condition, stop_int, tracing_terminated, & maxstep, istat) !! Performs integration along field line. If `condition` is present, !! performs field line tracing while `condition` is true real(FP), intent(in) :: x_init !! Start point (x-coordinate) real(FP), intent(in) :: y_init !! Start point (y-coordinate) real(FP), intent(in) :: phi_init !! Start point (phi-coordinate) real(FP), intent(in) :: dphi !! Toroidal distance to be traced class(equilibrium_t), target, intent(inout) :: equi !! Equilibrium defining the field line real(FP) , intent(out), optional :: x_end !! End point (x-coordinate) real(FP) , intent(out), optional :: y_end !! End point (y-coordinate) real(FP) , intent(out), optional :: phi_end !! Actual phi-coordinate end point of integration (relevant if !! `condition` is present) real(FP) , intent(out), optional :: arclen !! Length along field line real(FP) , intent(out), optional :: fluxexp !! Flux expansion, i.e. 1 / Btor integrated along fieldline procedure(condition_interface), optional :: condition !! Trace is performed until condition is true logical, intent(in), optional :: stop_int !! If = true, then the routine stops the integration once the !! condition is met. Otherwise, integration continues with integrand !! multiplied by zero. The default setting is false. logical, intent(out), optional :: tracing_terminated !! True, if given as an input, and the tracing was terminated due the !! fulfillment of the given condition. False by default. real(FP) , intent(in), optional :: maxstep !! Maximum stepsize for the integration integer, intent(out), optional :: istat !! PARALLAX error status of trace routine. If present, the trace routine !! will trigger a warning in case of error, and return an error code. If !! not present, the trace routine will trigger an error. type(dop853_equi_t) :: odeint logical :: odestat_ok integer :: idid, trace_error_code real(DP) :: phi_init_dp, phi_end_dp real(DP), dimension(neq) :: y if (present(istat)) then istat = PARALLAX_SUCCESS trace_error_code = PARALLAX_WRN_TRACE else trace_error_code = PARALLAX_ERR_TRACE endif if(present(tracing_terminated)) tracing_terminated = .false. ! For dphi = 0 if (abs(dphi) < FP_EPS) then if (present(x_end)) x_end = x_init if (present(y_end)) y_end = y_init if (present(phi_end)) phi_end = phi_init if (present(arclen)) arclen = 0.0_FP if (present(fluxexp)) fluxexp = 0.0_FP return endif ! Initialise ODE integrator call odeint%set_equilibrium(equi) if (present(condition)) then call odeint%set_condition(condition) endif if (present(maxstep)) then call odeint%initialize(n=neq, iprint=get_stderr(), hmax=maxstep, & status_ok=odestat_ok) else call odeint%initialize(n=neq, iprint=get_stderr(),& status_ok=odestat_ok) endif if (.not. odestat_ok) then call handle_error("DOP853 integrator initialization failed!", & trace_error_code, __LINE__, __FILE__) if (present(istat)) istat = PARALLAX_ERR_TRACE endif ! Set initial condition phi_init_dp = real(phi_init, kind=DP) phi_end_dp = real(phi_init + dphi, kind=DP) y(1) = phi_init_dp y(2) = real(x_init, kind=DP) y(3) = real(y_init, kind=DP) y(4) = 0.0_DP y(5) = 0.0_DP ! Perform integration if(present(stop_int)) then if(stop_int) then odeint%stop_integration = stop_int call odeint%integrate(phi_init_dp, y, phi_end_dp, rtol, atol, & iout=1, idid=idid) if (idid < 1) then call handle_error("Tracing failed, DOP853 did not converge!", & trace_error_code, __LINE__, __FILE__, & additional_info=& error_info_t("DOP853 status: ", [idid])) if (present(istat)) istat = PARALLAX_ERR_TRACE endif if (idid == 2) then if(present(tracing_terminated)) tracing_terminated = .true. endif else call odeint%integrate(phi_init_dp, y, phi_end_dp, rtol, atol, & iout=0, idid=idid) if (idid /= 1) then call handle_error("Tracing failed, DOP853 did not converge!", & trace_error_code, __LINE__, __FILE__, & additional_info=& error_info_t("DOP853 status: ", [idid])) if (present(istat)) istat = PARALLAX_ERR_TRACE endif endif else call odeint%integrate(phi_init_dp, y, phi_end_dp, rtol, atol, & iout=0, idid=idid) if (idid /= 1) then call handle_error("Tracing failed, DOP853 did not converge!", & trace_error_code, __LINE__, __FILE__, & additional_info=& error_info_t("DOP853 status: ", [idid])) if (present(istat)) istat = PARALLAX_ERR_TRACE endif endif ! Get solution if (present(phi_end)) phi_end = real(y(1), kind=FP) if (present(x_end)) x_end = real(y(2), kind=FP) if (present(y_end)) y_end = real(y(3), kind=FP) if (present(arclen)) arclen = real(abs(y(4)), kind=FP) if (present(fluxexp)) fluxexp = real(abs(y(5)), kind=FP) call odeint%destroy() end subroutine trace end module fieldline_tracer_m