diff --git a/src/beams.f90 b/src/beams.f90 index b6c8759..e85c25e 100644 --- a/src/beams.f90 +++ b/src/beams.f90 @@ -96,10 +96,10 @@ contains ! directions (mm⁻¹) ! - φ_w, φ_R are the angles of the amplitude and phase ellipses (deg) - use gray_params, only : antenna_parameters - use simplespline, only : spli, difcs - use utils, only : get_free_unit,locate - use logger, only : log_error + use gray_params, only : antenna_parameters + use splines, only : spline_simple + use utils, only : get_free_unit,locate + use logger, only : log_error implicit none @@ -108,13 +108,14 @@ contains integer, intent(in), optional :: unit ! local variables - integer :: u, iopt, ier, nisteer, i, k, ii + integer :: u, nisteer, i, k, ii real(wp_) :: steer, dal real(wp_), dimension(:), allocatable :: & alphastv, betastv, x00v, y00v, & - z00v, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v, & - cbeta, cx0, cy0, cz0, cwaist1, cwaist2, & - crci1, crci2, cphi1, cphi2 + z00v, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v + type(spline_simple) :: beta, waist1, waist2, & + rci1, rci2, phi1, phi2, & + x0, y0, z0 integer :: err u = get_free_unit(unit) @@ -128,13 +129,10 @@ contains read(u,*) params%fghz read(u,*) nisteer - allocate(alphastv(nisteer), betastv(nisteer), waist1v(nisteer), & - waist2v(nisteer), rci1v(nisteer), rci2v(nisteer), & - phi1v(nisteer), phi2v(nisteer), x00v(nisteer), & - y00v(nisteer), z00v(nisteer), cbeta(4*nisteer), & - cx0(4*nisteer), cy0(4*nisteer), cz0(4*nisteer), & - cwaist1(4*nisteer), cwaist2(4*nisteer), crci1(4*nisteer), & - crci2(4*nisteer), cphi1(4*nisteer), cphi2(4*nisteer)) + allocate(alphastv(nisteer), betastv(nisteer), waist1v(nisteer), & + waist2v(nisteer), rci1v(nisteer), rci2v(nisteer), & + phi1v(nisteer), phi2v(nisteer), x00v(nisteer), & + y00v(nisteer), z00v(nisteer)) do i=1,nisteer read(u, *) steer, alphastv(i), betastv(i), & @@ -151,34 +149,33 @@ contains z00v = 0.1_wp_ * z00v waist1v = 0.1_wp_ * waist1v waist2v = 0.1_wp_ * waist2v - rci1v = 10._wp_ * rci1v - rci2v = 10._wp_ * rci2v + rci1v = 10 * rci1v + rci2v = 10 * rci2v - iopt = 0 - call difcs(alphastv, betastv, nisteer, iopt, cbeta, ier) - call difcs(alphastv, waist1v, nisteer, iopt, cwaist1, ier) - call difcs(alphastv, rci1v, nisteer, iopt, crci1, ier) - call difcs(alphastv, waist2v, nisteer, iopt, cwaist2, ier) - call difcs(alphastv, rci2v, nisteer, iopt, crci2, ier) - call difcs(alphastv, phi1v, nisteer, iopt, cphi1, ier) - call difcs(alphastv, phi2v, nisteer, iopt, cphi2, ier) - call difcs(alphastv, x00v, nisteer, iopt, cx0, ier) - call difcs(alphastv, y00v, nisteer, iopt, cy0, ier) - call difcs(alphastv, z00v, nisteer, iopt, cz0, ier) + call beta%init(alphastv, betastv, nisteer) + call waist1%init(alphastv, waist1v, nisteer) + call waist2%init(alphastv, waist2v, nisteer) + call rci1%init(alphastv, rci1v, nisteer) + call rci2%init(alphastv, rci2v, nisteer) + call phi1%init(alphastv, phi1v, nisteer) + call phi2%init(alphastv, phi2v, nisteer) + call x0%init(alphastv, x00v, nisteer) + call y0%init(alphastv, y00v, nisteer) + call z0%init(alphastv, z00v, nisteer) if((params%alpha > alphastv(1)) .and. (params%alpha < alphastv(nisteer))) then call locate(alphastv, nisteer, params%alpha , k) dal = params%alpha - alphastv(k) - params%beta = spli(cbeta, nisteer, k, dal) - params%pos(1) = spli(cx0, nisteer, k, dal) - params%pos(2) = spli(cy0, nisteer, k, dal) - params%pos(3) = spli(cz0, nisteer, k, dal) - params%w(1) = spli(cwaist1, nisteer, k, dal) - params%w(2) = spli(cwaist2, nisteer, k, dal) - params%ri(1) = spli(crci1, nisteer, k, dal) - params%ri(2) = spli(crci2, nisteer, k, dal) - params%phi(1) = spli(cphi1, nisteer, k, dal) - params%phi(2) = spli(cphi2, nisteer, k, dal) + params%beta = beta%raw_eval(k, dal) + params%pos(1) = x0%raw_eval(k, dal) + params%pos(2) = y0%raw_eval(k, dal) + params%pos(3) = z0%raw_eval(k, dal) + params%w(1) = waist1%raw_eval(k, dal) + params%w(2) = waist2%raw_eval(k, dal) + params%ri(1) = rci1%raw_eval(k, dal) + params%ri(2) = rci2%raw_eval(k, dal) + params%phi(1) = phi1%raw_eval(k, dal) + params%phi(2) = phi2%raw_eval(k, dal) else ! params%alpha outside table range if(params%alpha >= alphastv(nisteer)) ii=nisteer @@ -196,10 +193,20 @@ contains params%phi(2) = phi2v(ii) end if - deallocate(alphastv, betastv, waist1v, waist2v, rci1v, rci2v, & - phi1v, phi2v, x00v, y00v, z00v, cbeta, & - cx0, cy0, cz0, cwaist1, cwaist2, & - crci1, crci2, cphi1, cphi2) + deallocate(alphastv, betastv, waist1v, waist2v, & + rci1v, rci2v, phi1v, phi2v, & + x00v, y00v, z00v) + + call beta%deinit + call waist1%deinit + call waist2%deinit + call rci1%deinit + call rci2%deinit + call phi1%deinit + call phi2%deinit + call x0%deinit + call y0%deinit + call z0%deinit end subroutine read_beam1 diff --git a/src/coreprofiles.f90 b/src/coreprofiles.f90 index 888b132..a9afd97 100644 --- a/src/coreprofiles.f90 +++ b/src/coreprofiles.f90 @@ -1,4 +1,4 @@ -! This modules handles the loading, interpolation and evaluation of the +! This module handles the loading, interpolation and evaluation of the ! plasma profiles (density, temperature, effective charge) ! ! Two kinds of profiles are supported: analytical (suffix `_an` in the @@ -6,20 +6,10 @@ ! the data is interpolated using splines. module coreprofiles use const_and_precisions, only : wp_, zero, one + use splines, only : spline_simple, spline_1d implicit none - ! Parameters of the plasma profiles splines - type spline_parameters - integer :: ndata ! Number of data points - integer :: nknots ! Number of spline knots - ! Density spline (ψ, knots, B-spline coefficients) - real(wp_), dimension(:), allocatable :: knots, coeffs - ! Temperature and effective charge arrays (ψ, T(ψ), Zeff(ψ)) - real(wp_), dimension(:), allocatable :: psi - real(wp_), dimension(:, :), allocatable :: temp, zeff - end type - ! Parameters of the C² polynomial tail of the density spline type density_tail real(wp_) :: start ! ψ₀, start of the tail @@ -38,10 +28,12 @@ module coreprofiles real(wp_) :: zeff ! Effective charge end type - ! Global variable storing the state of the module - type(spline_parameters), save :: spline - type(density_tail), save :: tail - type(analytic_model), save :: model + ! Global variables storing the state of the module + type(spline_1d), save :: dens_spline + type(spline_simple), save :: temp_spline + type(spline_simple), save :: zeff_spline + type(density_tail), save :: tail + type(analytic_model), save :: model private public read_profiles, read_profiles_an ! Reading data files @@ -58,7 +50,6 @@ contains ! ! Note: density has units of 10¹⁹ m⁻³. use gray_params, only : iprof - use dierckx, only : splev, splder use logger, only : log_error implicit none @@ -68,14 +59,11 @@ contains real(wp_), intent(out) :: dens, ddens ! density and first derivative ! local variables - integer :: ier ! dierck error code - real(wp_) :: f(1) ! dierck output (must be an array) - real(wp_) :: wrkfd(spline%ndata+4) ! dierck working space array - character(256) :: msg ! for log messages formatting + character(256) :: msg ! for log messages formatting ! Initialise both to zero - dens = zero - ddens = zero + dens = 0 + ddens = 0 ! Outside the tail end both density and its ! derivatives are identically zero @@ -97,16 +85,10 @@ contains ! Use the interpolating spline when in range ! Evaluate the spline - ier = 0 - call splev(spline%knots, spline%nknots, spline%coeffs, & - 3, [psin], f, 1, ier) - dens = f(1) + dens = dens_spline%eval(psin) + ddens = dens_spline%deriv(psin) ! Evaluate the spline 1st derivative - ier = 0 - call splder(spline%knots, spline%nknots, spline%coeffs, & - 3, 1, [psin], f, 1, wrkfd, ier) - ddens = f(1) if (abs(dens) < 1.0e-10_wp_) dens = zero else ! Use a C² polynomial extension outside (ψ > ψ₀) @@ -158,9 +140,7 @@ contains ! normalised poloidal flux. ! ! Note: temperature has units of keV. - use gray_params, only : iprof - use utils, only : locate - use simplespline, only : spli + use gray_params, only : iprof implicit none @@ -169,10 +149,9 @@ contains real(wp_) :: temp ! local variables - integer :: k - real(wp_) :: proft, dps + real(wp_) :: proft - temp = zero + temp = 0 if (psin >= 1 .or. psin < 0) return if (iprof == 0) then ! Use the analytical model @@ -183,10 +162,7 @@ contains temp = (model%te0 - model%te1)*proft + model%te1 else ! Use the interpolated numerical data - call locate(spline%psi, spline%ndata, psin, k) - k = max(1, min(k, spline%ndata - 1)) - dps = psin - spline%psi(k) - temp = spli(spline%temp, spline%ndata, k, dps) + temp = temp_spline%eval(psin) endif end function temp @@ -195,8 +171,6 @@ contains ! Computes the effective charge Zeff as a ! function of the normalised poloidal flux. use gray_params, only : iprof - use utils, only : locate - use simplespline, only : spli implicit none @@ -204,10 +178,6 @@ contains real(wp_), intent(in) :: psin real(wp_) :: fzeff - ! local variables - integer :: k - real(wp_) :: dps - fzeff = one if (psin >= 1 .or. psin < 0) return if (iprof == 0) then @@ -215,10 +185,7 @@ contains fzeff = model%zeff else ! Use the interpolated numerical data - call locate(spline%psi, spline%ndata, psin, k) - k = max(1, min(k, spline%ndata - 1)) - dps = psin - spline%psi(k) - fzeff = spli(spline%zeff, spline%ndata, k, dps) + fzeff = zeff_spline%eval(psin) endif end function fzeff @@ -378,8 +345,6 @@ contains ! When `launch_pos` (cartesian launch coordinates in cm) is present, ! the subroutine will also check that the wave launcher is strictly ! outside the reconstructed plasma density boundary. - use simplespline, only : difcs - use dierckx, only : curfit, splev, splder use gray_params, only : profiles_parameters, profiles_data use logger, only : log_debug, log_info, log_warning, log_error @@ -390,68 +355,28 @@ contains type(profiles_data), intent(inout) :: data real(wp_), optional, intent(in) :: launch_pos(3) - ! curfit parameters - integer, parameter :: iopt = 0 ! smoothing spline mode - integer, parameter :: kspl = 3 ! order of spline (cubic) - ! local variables - integer :: n, npest, ier - real(wp_) :: xb, xe, fp, ssplne_loc - - ! working space arrays for the dierckx functions - integer :: lwrkf - real(wp_), dimension(:), allocatable :: wf, wrkf - integer, dimension(:), allocatable :: iwrkf + integer :: n, err ! for log messages formatting character(256) :: msg n = size(data%psrad) - npest = n + 4 - lwrkf = n*4 + npest*16 - allocate(wrkf(lwrkf), iwrkf(npest), wf(n)) - - ssplne_loc=params%sspld - - ! If necessary, reallocate the spline arrays - if (.not. allocated(spline%psi)) then - allocate(spline%psi(n), spline%temp(n, 4), spline%zeff(n, 4)) - else - if (size(spline%psi) < n) then - deallocate(spline%psi, spline%temp, spline%zeff) - allocate(spline%psi(n), spline%temp(n, 4), spline%zeff(n, 4)) - end if - end if - if (.not. allocated(spline%coeffs)) then - allocate(spline%knots(npest), spline%coeffs(npest)) - else - if (size(spline%coeffs) < npest) then - deallocate(spline%knots, spline%coeffs) - allocate(spline%knots(npest), spline%coeffs(npest)) - end if - end if ! Spline interpolation of temperature and effective charge - call difcs(data%psrad, data%terad, n, iopt, spline%temp, ier) - call difcs(data%psrad, data%zfc, n, iopt, spline%zeff, ier) - spline%psi = data%psrad - spline%ndata = n + call temp_spline%init(data%psrad, data%terad, n) + call zeff_spline%init(data%psrad, data%zfc, n) - ! Spline interpolation of density - xb = zero - xe = data%psrad(n) - wf(:) = one - call curfit(iopt, n, data%psrad, data%derad, wf, xb, xe, kspl, & - ssplne_loc, npest, spline%nknots, spline%knots, & - spline%coeffs, fp, wrkf, lwrkf, iwrkf, ier) - ! if ier=-1 data are re-fitted using sspl=0 - if (ier == -1) then - call log_warning('curfit failed with error -1: re-fitting with '// & - 's=0', mod='coreprofiles', proc='density') - ssplne_loc = zero - call curfit(iopt, n, data%psrad, data%derad, wf, xb, xe, kspl, & - ssplne_loc, npest, spline%nknots, spline%knots, & - spline%coeffs, fp, wrkf, lwrkf, iwrkf, ier) + ! Spline interpolation of density (smooth spline) + call dens_spline%init(data%psrad, data%derad, n, range=[zero, data%psrad(n)], & + tension=params%sspld, err=err) + + ! if failed, re-fit with an interpolating spline (zero tension) + if (err == -1) then + call log_warning('density fit failed with error -1: re-fitting with '// & + 'zero tension', mod='coreprofiles', proc='density') + call dens_spline%init(data%psrad, data%derad, n, & + range=[zero, data%psrad(n)], tension=zero) end if ! Computation of the polynomial tail parameters @@ -460,9 +385,9 @@ contains ! at the edge. The spline thus has to be extended to transition ! smoothly from the last profile point to 0 outside the plasma. block - real(wp_), dimension(1) :: s0, s1, s2 ! spline, 1st, 2nd derivative - real(wp_), dimension(1) :: delta4 ! discriminant Δ/4 of q(x) - real(wp_), dimension(1) :: x0, x1 ! vertex of q(x), solution + real(wp_) :: s0, s1, s2 ! spline, 1st, 2nd derivative + real(wp_) :: delta4 ! discriminant Δ/4 of q(x) + real(wp_) :: x0, x1 ! vertex of q(x), solution ! Compute the coefficients of a 2nd order Taylor polinomial to ! extend the spline beyond the last point: @@ -471,12 +396,9 @@ contains ! ! where s(ψ) is the spline and ψ₀ the last point. ! - call splev(spline%knots, spline%nknots, spline%coeffs, kspl, & - data%psrad(n:n), s0, 1, ier) - call splder(spline%knots, spline%nknots, spline%coeffs, kspl, 1, & - data%psrad(n:n), s1, 1, wrkf(1:spline%nknots), ier) - call splder(spline%knots, spline%nknots, spline%coeffs, kspl, 2, & - data%psrad(n:n), s2, 1, wrkf(1:spline%nknots), ier) + s0 = dens_spline%eval(data%psrad(n)) + s1 = dens_spline%deriv(data%psrad(n), order=1) + s2 = dens_spline%deriv(data%psrad(n), order=2) ! Determine where to end the tail (to ensure the density remains ! positive) from the zeros of the Taylor polynomial p(ψ) @@ -491,7 +413,7 @@ contains x0 = -s1 / s2 ! vertex of parabola y=q(x) delta4 = (s1 / s2)**2 - 2*s0/s2 ! Δ/4 of q(x) - if (delta4(1) > 0) then + if (delta4 > 0) then ! Pick the smallest positive zero (implying >ψ₀) x1 = x0 + sign(sqrt(delta4), sqrt(delta4) - x0) else @@ -503,10 +425,10 @@ contains ! Store the tail parameters tail%start = data%psrad(n) - tail%end = tail%start + x1(1) - tail%value = s0(1) - tail%deriv1 = s1(1) - tail%deriv2 = s2(1) + tail%end = tail%start + x1 + tail%value = s0 + tail%deriv1 = s1 + tail%deriv2 = s2 end block ! Make sure the wave launcher does not fall inside the tail @@ -550,8 +472,6 @@ contains write (msg, '(a,g0.4)') 'density boundary: ψ=', tail%end call log_info(msg, mod='coreprofiles', proc='set_profiles_spline') - - deallocate(iwrkf, wrkf, wf) end subroutine set_profiles_spline @@ -560,11 +480,9 @@ contains implicit none - if (allocated(spline%psi)) deallocate(spline%psi) - if (allocated(spline%temp)) deallocate(spline%temp) - if (allocated(spline%zeff)) deallocate(spline%zeff) - if (allocated(spline%knots)) deallocate(spline%knots) - if (allocated(spline%coeffs)) deallocate(spline%coeffs) + call dens_spline%deinit + call temp_spline%deinit + call zeff_spline%deinit end subroutine unset_profiles_spline diff --git a/src/dierckx.f90 b/src/dierckx.f90 index 6d1f484..54a349b 100644 --- a/src/dierckx.f90 +++ b/src/dierckx.f90 @@ -4606,4 +4606,4 @@ contains end do end subroutine fpcuro -end module dierckx \ No newline at end of file +end module dierckx diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index 3bae324..7feb2b9 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -1,36 +1,67 @@ +! This module handles the loading, interpolation and evaluation of the +! MHD equilibrium data (poloidal current function, safety factor, +! poloidal flux, magnetic field, etc.) +! +! Two kinds of data are supported: analytical (suffix `_an` in the +! subroutine names) or numerical (suffix `_spline`). For the latter, the +! the data is interpolated using splines. module equilibrium use const_and_precisions, only : wp_ + use splines, only : spline_simple, spline_1d, spline_2d + implicit none - - real(wp_), save :: btaxis,rmaxis,zmaxis,sgnbphi - real(wp_), save :: btrcen ! used only for jcd_astra def. - real(wp_), save :: rcen ! computed as fpol(a)/btrcen - real(wp_), save :: rmnm,rmxm,zmnm,zmxm - real(wp_), save :: zbinf,zbsup -! real(wp_), save :: rup,zup,rlw,zlw - integer, parameter :: kspl=3,ksplp=kspl+1 + ! Parameters of the analytical equilibrium model + type analytic_model + real(wp_) :: q0 ! Safety factor at the magnetic axis + real(wp_) :: q1 ! Safety factor at the edge + real(wp_) :: lq ! Exponent for the q(ρ) power law + end type -! === 2d spline psi(r,z), normalization and derivatives ========== - integer, save :: nsr, nsz - real(wp_), save :: psia, psiant, psinop, psnbnd - real(wp_), dimension(:), allocatable, save :: tr,tz - real(wp_), dimension(:), allocatable, save :: cceq, cceq01, cceq10, & - cceq20, cceq02, cceq11 + ! Order of the splines + integer, parameter :: kspl=3, ksplp=kspl + 1 -! === 1d spline fpol(psi) ======================================== -! integer, save :: npsiest - integer, save :: nsf - real(wp_), dimension(:), allocatable, save :: tfp, cfp + ! Global variable storing the state of the module + + ! Splines + type(spline_1d), save :: fpol_spline ! Poloidal current function F(ψ) + type(spline_2d), save :: psi_spline ! Poloidal flux ψ(R,z) + type(spline_simple), save :: q_spline ! Safey factor q(ψ) + type(spline_simple), save :: rhop_spline, rhot_spline ! Normalised radii ρ_p(ρ_t), ρ_t(ρ_p) + + ! Analytic model + type(analytic_model), save :: model + + ! More parameters real(wp_), save :: fpolas + real(wp_), save :: psia, psiant, psinop + real(wp_), save :: phitedge, aminor + real(wp_), save :: btaxis,rmaxis,zmaxis,sgnbphi + real(wp_), save :: btrcen ! used only for jcd_astra def. + real(wp_), save :: rcen ! computed as fpol(a)/btrcen + real(wp_), save :: rmnm, rmxm, zmnm, zmxm + real(wp_), save :: zbinf, zbsup -! === 1d spline rhot(rhop), rhop(rhot), q(psi) =================== -! computed on psinr,rhopnr [,rhotnr] arrays - integer, save :: nq,nrho - real(wp_), dimension(:), allocatable, save :: psinr,rhopr,rhotr - real(wp_), dimension(:,:), allocatable, save :: cq,crhop,crhot - real(wp_), save :: phitedge,aminor - real(wp_), save :: q0,qa,alq + private + public read_eqdsk, read_equil_an ! Reading data files + public scale_equil, change_cocos ! Transforming data + public equian ! Analytical model + public equinum_psi, equinum_fpol ! Interpolated data + public fq, bfield, tor_curr, tor_curr_psi ! Accessing local quantities + public frhopol, frhopolv, frhotor ! Converting between poloidal/toroidal flux + public set_equil_spline, set_equil_an ! Initialising internal state + public unset_equil_spline ! Deinitialising internal state + + ! Members exposed for magsurf_data + public kspl, psi_spline, q_spline, points_tgo + + ! Members exposed to gray_core and more + public psia, psiant, psinop + public phitedge, aminor + public btaxis,rmaxis,zmaxis,sgnbphi + public btrcen, rcen + public rmnm, rmxm, zmnm, zmxm + public zbinf, zbsup contains @@ -181,7 +212,7 @@ contains ! local variables integer :: i, u, nlim integer :: err - real(wp_) :: rr0m, zr0m, rpam, b0, q0, qa, alq + real(wp_) :: rr0m, zr0m, rpam, b0 u = get_free_unit(unit) @@ -194,7 +225,7 @@ contains read(u, *) rr0m, zr0m, rpam read(u, *) b0 - read(u, *) q0, qa, alq + read(u, *) model%q0, model%q1, model%lq if(allocated(data%rv)) deallocate(data%rv) if(allocated(data%zv)) deallocate(data%zv) @@ -205,7 +236,7 @@ contains data%rv = [rr0m, rpam] data%zv = [zr0m] data%fpol = [b0*rr0m] - data%qpsi = [q0, qa, alq] + data%qpsi = [model%q0, model%q1, model%lq] if(ipass >= 2) then ! get size of boundary and limiter arrays and allocate them @@ -278,6 +309,7 @@ contains subroutine decode_cocos(cocos, exp2pi, phiccw, psiincr, qpos) + ! Extracts the sign and units conventions from a COCOS index implicit none ! subroutine arguments @@ -305,7 +337,7 @@ contains end subroutine decode_cocos - subroutine eq_scal(params, data) + subroutine scale_equil(params, data) ! Rescale the magnetic field (B) and the plasma current (I) ! and/or force their signs. ! @@ -339,7 +371,7 @@ contains data%fpol = data%fpol * params%factb params%sgni = int(sign(one, -data%psia)) params%sgnb = int(sign(one, last_fpol)) - end subroutine eq_scal + end subroutine scale_equil subroutine set_equil_spline(params, data) @@ -347,7 +379,6 @@ contains ! in their respective global variables, see the top of this file. use const_and_precisions, only : zero, one use gray_params, only : equilibrium_parameters, equilibrium_data - use dierckx, only : regrid, coeff_parder, curfit, splev use gray_params, only : iequil use reflections, only : inside use utils, only : vmaxmin, vmaxmini @@ -359,55 +390,44 @@ contains type(equilibrium_parameters), intent(in) :: params type(equilibrium_data), intent(in) :: data - ! local constants - integer, parameter :: iopt=0 - ! local variables - integer :: liwrk,lwrk,lw10,lw01,lw20,lw02,lw11,lwrkf - integer :: nr,nz,nrest,nzest,npsest,nrz,npsi,nbnd,ibinf,ibsup - real(wp_) :: sspln,fp,rax0,zax0,psinoptmp,psinxptmp - real(wp_) :: rbmin,rbmax,rbinf,rbsup,r1,z1 - real(wp_), dimension(size(data%psinr)) :: rhotn - real(wp_), dimension(1) :: fpoli - real(wp_), dimension(:), allocatable :: rv1d,zv1d,fvpsi,wf,wrk - integer, dimension(:), allocatable :: iwrk - integer :: ier,ixploc,info,i,j,ij + integer :: nr, nz, nrest, nzest, npsest, nrz, npsi, nbnd, ibinf, ibsup + real(wp_) :: tension, rax0, zax0, psinoptmp, psinxptmp + real(wp_) :: rbmin, rbmax, rbinf, rbsup, r1, z1 + real(wp_), dimension(size(data%psinr)) :: rhotn + real(wp_), dimension(:), allocatable :: rv1d, zv1d, fvpsi, wf + integer :: ier, ixploc, info, i, j, ij character(256) :: msg ! for log messages formatting - ! compute array sizes and prepare working space arrays - nr=size(data%rv) - nz=size(data%zv) - nrz=nr*nz - nrest=nr+ksplp - nzest=nz+ksplp - lwrk=4+nrest*nz+(nrest+nzest)*(2*kspl+5)+(nr+nz)*ksplp+max(nz,nrest) - liwrk=nz+nr+nrest+nzest+kspl - - npsi=size(data%psinr) - npsest=npsi+ksplp - lwrkf=npsi*ksplp+npsest*(7+3*kspl) - - allocate(wrk(max(lwrk,lwrkf)),iwrk(max(liwrk,npsest))) - - ! spline fitting/interpolation of data%psin(i,j) and derivatives - - ! allocate knots and spline coefficients arrays - if (allocated(tr)) deallocate(tr) - if (allocated(tz)) deallocate(tz) - if (allocated(cceq)) deallocate(cceq) - allocate(tr(nrest),tz(nzest),cceq(nrz)) + ! compute array sizes + nr = size(data%rv) + nz = size(data%zv) + nrz = nr*nz + npsi = size(data%psinr) + nrest = nr + ksplp + nzest = nz + ksplp + npsest = npsi + ksplp ! length in m !!! - rmnm=data%rv(1) - rmxm=data%rv(nr) - zmnm=data%zv(1) - zmxm=data%zv(nz) + rmnm = data%rv(1) + rmxm = data%rv(nr) + zmnm = data%zv(1) + zmxm = data%zv(nz) + + ! Spline interpolation of ψ(R, z) if (iequil>2) then ! data valid only inside boundary (data%psin=0 outside), e.g. source==ESCO ! presence of boundary anticipated here to filter invalid data nbnd = min(size(data%rbnd), size(data%zbnd)) + ! allocate knots and spline coefficients arrays + if (allocated(psi_spline%knots_x)) deallocate(psi_spline%knots_x) + if (allocated(psi_spline%knots_y)) deallocate(psi_spline%knots_y) + if (allocated(psi_spline%coeffs)) deallocate(psi_spline%coeffs) + allocate(psi_spline%knots_x(nrest), psi_spline%knots_y(nzest)) + allocate(psi_spline%coeffs(nrz)) + ! determine number of valid grid points nrz=0 do j=1,nz @@ -439,86 +459,79 @@ contains end do end do - ! fit as a scattered set of points + ! Fit as a scattered set of points ! use reduced number of knots to limit memory comsumption ? - nsr=nr/4+4 - nsz=nz/4+4 - sspln=params%ssplps - call scatterspl(rv1d,zv1d,fvpsi,wf,nrz,kspl,sspln, & - rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cceq,ier) - ! if ier=-1 data are fitted using params%ssplps=0 - if(ier.eq.-1) then - sspln=0.0_wp_ - nsr=nr/4+4 - nsz=nz/4+4 - call scatterspl(rv1d,zv1d,fvpsi,wf,nrz,kspl,sspln, & - rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cceq,ier) + psi_spline%nknots_x=nr/4+4 + psi_spline%nknots_y=nz/4+4 + tension = params%ssplps + call scatterspl(rv1d, zv1d, fvpsi, wf, nrz, kspl, tension, & + rmnm, rmxm, zmnm, zmxm, & + psi_spline%knots_x, psi_spline%nknots_x, & + psi_spline%knots_y, psi_spline%nknots_y, & + psi_spline%coeffs, ier) + + ! if failed, re-fit with an interpolating spline (zero tension) + if(ier == -1) then + tension = 0 + psi_spline%nknots_x=nr/4+4 + psi_spline%nknots_y=nz/4+4 + call scatterspl(rv1d, zv1d, fvpsi, wf, nrz, kspl, tension, & + rmnm, rmxm, zmnm, zmxm, & + psi_spline%knots_x, psi_spline%nknots_x, & + psi_spline%knots_y, psi_spline%nknots_y, & + psi_spline%coeffs, ier) end if - deallocate(rv1d,zv1d,wf,fvpsi) + deallocate(rv1d, zv1d, wf, fvpsi) ! reset nrz to the total number of grid points for next allocations - nrz=nr*nz + nrz = nr*nz else ! iequil==2: data are valid on the full R,z grid - ! reshape 2D psi array to 1D (transposed) array and compute spline coeffs + ! reshape 2D ψ array to 1D (transposed) allocate(fvpsi(nrz)) - fvpsi=reshape(transpose(data%psin),(/nrz/)) - sspln=params%ssplps - call regrid(iopt,nr,data%rv,nz,data%zv,fvpsi,rmnm,rmxm,zmnm,zmxm, & - kspl,kspl,sspln,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, & - wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier) - ! if ier=-1 data are re-fitted using params%ssplps=0 - if(ier==-1) then - sspln=0.0_wp_ - call regrid(iopt,nr,data%rv,nz,data%zv,fvpsi,rmnm,rmxm,zmnm,zmxm, & - kspl,kspl,sspln,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, & - wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier) + fvpsi = reshape(transpose(data%psin), [nrz]) + + ! compute spline coefficients + call psi_spline%init(data%rv, data%zv, fvpsi, nr, nz, & + range=[rmnm, rmxm, zmnm, zmxm], & + tension=params%ssplps, err=ier) + + ! if failed, re-fit with an interpolating spline (zero tension) + if(ier == -1) then + call psi_spline%init(data%rv, data%zv, fvpsi, nr, nz, & + range=[rmnm, rmxm, zmnm, zmxm], & + tension=zero) end if deallocate(fvpsi) end if - ! compute spline coefficients for psi partial derivatives - lw10 = nr*(ksplp-1) + nz*ksplp + nrz - lw01 = nr*ksplp + nz*(ksplp-1) + nrz - lw20 = nr*(ksplp-2) + nz*ksplp + nrz - lw02 = nr*ksplp + nz*(ksplp-2) + nrz - lw11 = nr*(ksplp-1) + nz*(ksplp-1) + nrz - if (allocated(cceq10)) deallocate(cceq10) - if (allocated(cceq01)) deallocate(cceq01) - if (allocated(cceq20)) deallocate(cceq20) - if (allocated(cceq02)) deallocate(cceq02) - if (allocated(cceq11)) deallocate(cceq11) - allocate(cceq10(lw10),cceq01(lw01),cceq20(lw20),cceq02(lw02),cceq11(lw11)) - call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,1,0,cceq10,lw10,ier) - call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,0,1,cceq01,lw01,ier) - call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,2,0,cceq20,lw20,ier) - call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,0,2,cceq02,lw02,ier) - call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,1,1,cceq11,lw11,ier) + ! compute spline coefficients for ψ(R,z) partial derivatives + call psi_spline%init_deriv(nr, nz, 1, 0) ! ∂ψ/∂R + call psi_spline%init_deriv(nr, nz, 0, 1) ! ∂ψ/∂z + call psi_spline%init_deriv(nr, nz, 1, 1) ! ∂²ψ/∂R∂z + call psi_spline%init_deriv(nr, nz, 2, 0) ! ∂²ψ/∂R² + call psi_spline%init_deriv(nr, nz, 0, 2) ! ∂²ψ/∂z² - ! Spline interpolation of data%fpol(i) + ! Spline interpolation of F(ψ) - ! Allocate knots and spline coefficients arrays - if (allocated(tfp)) deallocate(tfp) - if (allocated(cfp)) deallocate(cfp) - allocate(tfp(npsest),cfp(npsest)) + ! give a small weight to the last point allocate(wf(npsi)) - wf(1:npsi-1)=one - wf(npsi)=1.0e2_wp_ - call curfit(iopt,npsi,data%psinr,data%fpol,wf,zero,one,kspl,params%ssplf,npsest,nsf, & - tfp,cfp,fp,wrk(1:lwrkf),lwrkf,iwrk(1:npsest),ier) - call splev(tfp,nsf,cfp,kspl,data%psinr(npsi:npsi),fpoli,1,ier) - ! Set vacuum value used outside 0<=data%psin<=1 range - fpolas=fpoli(1) - sgnbphi=sign(one,fpolas) - ! Free temporary arrays - deallocate(wrk,iwrk,wf) + wf(1:npsi-1) = 1 + wf(npsi) = 1.0e2_wp_ + call fpol_spline%init(data%psinr, data%fpol, npsi, range=[zero, one], & + weights=wf, tension=params%ssplf) + deallocate(wf) - ! Re-normalize psi after spline computation + ! set vacuum value used outside 0 ≤ ψ ≤ 1 range + fpolas = fpol_spline%eval(data%psinr(npsi)) + sgnbphi = sign(one ,fpolas) + + ! Re-normalize ψ after spline computation ! Start with un-corrected psi - psia=data%psia - psinop=0.0_wp_ - psiant=1.0_wp_ + psia = data%psia + psinop = 0 + psiant = 1 ! Use provided boundary to set an initial guess ! for the search of O/X points @@ -604,7 +617,7 @@ contains ! Save Bt value on axis (required in flux_average and used in Jcd def) ! and vacuum value B0 at ref. radius data%rvac (used in Jcd_astra def) - call equinum_fpol(0.0_wp_,btaxis) + call equinum_fpol(zero, btaxis) btaxis = btaxis/rmaxis btrcen = fpolas/data%rvac rcen = data%rvac @@ -612,19 +625,23 @@ contains write (msg, '(2(a,g0.3))') 'Bt_center=', btrcen, ' Bt_axis=', btaxis call log_info(msg, mod='equilibrium', proc='set_equil_spline') - ! Compute rho_pol/rho_tor mapping based on input q profile - call setqphi_num(data%psinr,abs(data%qpsi),abs(psia),rhotn) - call set_rhospl(sqrt(data%psinr),rhotn) + ! Compute ρ_p/ρ_t mapping based on the input q profile + call setqphi_num(data%psinr, abs(data%qpsi), abs(psia), rhotn) + call set_rho_spline(sqrt(data%psinr), rhotn) end subroutine set_equil_spline subroutine scatterspl(x,y,z,w,n,kspl,sspl,xmin,xmax,ymin,ymax, & tx,nknt_x,ty,nknt_y,coeff,ierr) - use const_and_precisions, only : wp_, comp_eps - use dierckx, only : surfit + ! Computes the spline interpolation of a surface when + ! the data points are irregular, i.e. not on a uniform grid + use const_and_precisions, only : comp_eps + use dierckx, only : surfit + implicit none -! arguments + + ! subroutine arguments integer, intent(in) :: n real(wp_), dimension(n), intent(in) :: x, y, z real(wp_), dimension(n), intent(in) :: w @@ -636,7 +653,8 @@ contains integer, intent(inout) :: nknt_x, nknt_y real(wp_), dimension(nknt_x*nknt_y), intent(out) :: coeff integer, intent(out) :: ierr -! local variables + + ! local variables integer :: iopt real(wp_) :: resid integer :: u,v,km,ne,b1,b2,lwrk1,lwrk2,kwrk,nxest,nyest @@ -667,38 +685,33 @@ contains end subroutine scatterspl - subroutine setqphi_num(psinq,q,psia,rhotn) + ! Computes the spline of the safety factor q(ψ) use const_and_precisions, only : pi - use simplespline, only : difcs + implicit none -! arguments + + ! subroutine arguments real(wp_), dimension(:), intent(in) :: psinq,q real(wp_), intent(in) :: psia real(wp_), dimension(:), intent(out), optional :: rhotn -! local variables + + ! local variables real(wp_), dimension(size(q)) :: phit real(wp_) :: dx - integer, parameter :: iopt=0 - integer :: k,ier + integer :: k - nq=size(q) - if(allocated(psinr)) deallocate(psinr) - if(allocated(cq)) deallocate(cq) - allocate(psinr(nq),cq(nq,4)) + call q_spline%init(psinq, q, size(q)) - psinr=psinq - call difcs(psinr,q,nq,iopt,cq,ier) - -! Toroidal flux phi = 2*pi*Integral q dpsi - phit(1)=0.0_wp_ - do k=1,nq-1 - dx=psinr(k+1)-psinr(k) - phit(k+1)=phit(k) + dx*(cq(k,1) + dx*(cq(k,2)/2.0_wp_ + & - dx*(cq(k,3)/3.0_wp_ + dx* cq(k,4)/4.0_wp_) ) ) + ! Toroidal flux φ = 2π ∫q(ψ)dψ + phit(1)=0 + do k=1,q_spline%ndata-1 + dx=q_spline%data(k+1)-q_spline%data(k) + phit(k+1)=phit(k) + dx*(q_spline%coeffs(k,1) + dx*(q_spline%coeffs(k,2)/2 + & + dx*(q_spline%coeffs(k,3)/3 + dx* q_spline%coeffs(k,4)/4) ) ) end do - phitedge=phit(nq) - if(present(rhotn)) rhotn(1:nq)=sqrt(phit/phitedge) + phitedge=phit(q_spline%ndata) + if(present(rhotn)) rhotn(1:q_spline%ndata)=sqrt(phit/phitedge) phitedge=2*pi*psia*phitedge end subroutine setqphi_num @@ -730,348 +743,283 @@ contains q1 = data%qpsi(2) qexp = data%qpsi(3) - btaxis=bax - rmaxis=rax - zmaxis=zax - btrcen=bax - rcen=rax - aminor=a - zbinf=zmaxis-a - zbsup=zmaxis+a - q0=qax - qa=q1 - alq=qexp - sgnbphi=sign(one,bax) + btaxis = bax + rmaxis = rax + zmaxis = zax + btrcen = bax + rcen = rax + aminor = a + zbinf = zmaxis-a + zbsup = zmaxis+a + model%q0 = qax + model%q1 = q1 + model%lq = qexp + sgnbphi = sign(one, bax) - rmxm=rmaxis+aminor - rmnm=rmaxis-aminor - zmxm=zbsup - zmnm=zbinf + rmxm = rmaxis+aminor + rmnm = rmaxis-aminor + zmxm = zbsup + zmnm = zbinf if (present(n)) then - nq=n + q_spline%ndata = n else - nq=nqdef + q_spline%ndata = nqdef end if - if (allocated(psinr)) deallocate(psinr) - allocate(psinr(nq),rhotn(nq),rhopn(nq)) - - dr=one/(nq-1) - rhotn(1)=zero - psinr(1)=zero - res=zero - fq0=zero - do i=2,nq - rn=(i-1)*dr - qq=q0+(q1-q0)*rn**qexp - fq1=rn/qq - res=res+0.5_wp_*(fq1+fq0)*dr - fq0=fq1 - rhotn(i)=rn - psinr(i)=res + + if (allocated(q_spline%data)) deallocate(q_spline%data) + allocate(q_spline%data(q_spline%ndata)) + allocate(rhotn(q_spline%ndata), rhopn(q_spline%ndata)) + + dr = one/(q_spline%ndata - 1) + rhotn(1) = zero + q_spline%data(1) = zero + res = zero + fq0 = zero + do i = 2, q_spline%ndata + rn = (i - 1)*dr + qq = model%q0 + (model%q1 - model%q0) * rn**model%lq + fq1 = rn / qq + res = res + (fq1 + fq0)/2 * dr + fq0 = fq1 + rhotn(i) = rn + q_spline%data(i) = res end do - phitedge=btaxis*aminor**2 ! temporary - psia=res*phitedge - phitedge=pi*phitedge ! final - psinr=psinr/res - rhopn=sqrt(psinr) - - call set_rhospl(rhopn,rhotn) + phitedge = btaxis*aminor**2 ! temporary + psia = res*phitedge + phitedge = pi*phitedge ! final + q_spline%data = q_spline%data/res + rhopn = sqrt(q_spline%data) + + call set_rho_spline(rhopn, rhotn) end subroutine set_equil_an - subroutine set_rhospl(rhop,rhot) - use simplespline, only : difcs + subroutine set_rho_spline(rhop, rhot) + ! Computes the splines for converting between the poloidal (ρ_p) + ! and toroidal (ρ_t) normalised radii + implicit none -! arguments + + ! subroutine arguments real(wp_), dimension(:), intent(in) :: rhop, rhot -! local variables - integer, parameter :: iopt=0 - integer :: ier - - nrho=size(rhop) - if(allocated(rhopr)) deallocate(rhopr) - if(allocated(rhotr)) deallocate(rhotr) - if(allocated(crhop)) deallocate(crhop) - if(allocated(crhot)) deallocate(crhot) - allocate(rhopr(nrho),rhotr(nrho),crhop(nrho,4),crhot(nrho,4)) - - rhopr=rhop - rhotr=rhot - call difcs(rhotr,rhopr,nrho,iopt,crhop,ier) - call difcs(rhopr,rhotr,nrho,iopt,crhot,ier) - end subroutine set_rhospl + call rhop_spline%init(rhot, rhop, size(rhop)) + call rhot_spline%init(rhop, rhot, size(rhot)) + end subroutine set_rho_spline + subroutine equinum_psi(R, z, psi, dpsidr, dpsidz, & + ddpsidrr, ddpsidzz, ddpsidrz) + ! Computes the normalised poloidal flux ψ (and its derivatives) + ! as a function of (R, z) using the numerical equilibrium - subroutine equinum_psi(rpsim,zpsim,psinv,dpsidr,dpsidz, & - ddpsidrr,ddpsidzz,ddpsidrz) - use dierckx, only : fpbisp implicit none -! local constants - integer, parameter :: lwrk=8,liwrk=2 -! arguments - real(wp_), intent(in) :: rpsim,zpsim - real(wp_), intent(out), optional :: psinv,dpsidr,dpsidz, & - ddpsidrr,ddpsidzz,ddpsidrz -! local variables - integer, dimension(liwrk) :: iwrk - real(wp_), dimension(1) :: rrs,zzs,ffspl - real(wp_), dimension(lwrk) :: wrk + ! subroutine arguments + real(wp_), intent(in) :: R, z + real(wp_), intent(out), optional :: psi, dpsidr, dpsidz + real(wp_), intent(out), optional :: ddpsidrr, ddpsidzz, ddpsidrz -! here lengths are measured in meters + ! Note: here lengths are measured in meters + if (R <= rmxm .and. R >= rmnm .and. & + z <= zmxm .and. z >= zmnm) then - if (rpsim.le.rmxm .and. rpsim.ge.rmnm .and. & - zpsim.le.zmxm .and. zpsim.ge.zmnm) then - - if (present(psinv)) then - rrs(1)=rpsim - zzs(1)=zpsim - call fpbisp(tr,nsr,tz,nsz,cceq,3,3,rrs,1,zzs,1,ffspl, & - wrk(1),wrk(5),iwrk(1),iwrk(2)) - psinv=(ffspl(1)-psinop)/psiant - end if - if (present(dpsidr)) then - call sub_derpsi(rpsim,zpsim,1,0,dpsidr,cceq10) - end if - if (present(dpsidz)) then - call sub_derpsi(rpsim,zpsim,0,1,dpsidz,cceq01) - end if - if (present(ddpsidrr)) then - call sub_derpsi(rpsim,zpsim,2,0,ddpsidrr,cceq20) - end if - if (present(ddpsidzz)) then - call sub_derpsi(rpsim,zpsim,0,2,ddpsidzz,cceq02) - end if - if (present(ddpsidrz)) then - call sub_derpsi(rpsim,zpsim,1,1,ddpsidrz,cceq11) - end if + if (present(psi)) psi = (psi_spline%eval(R, z) - psinop) / psiant + if (present(dpsidr)) dpsidr = psi_spline%deriv(R, z, 1, 0) * psia + if (present(dpsidz)) dpsidz = psi_spline%deriv(R, z, 0, 1) * psia + if (present(ddpsidrr)) ddpsidrr = psi_spline%deriv(R, z, 2, 0) * psia + if (present(ddpsidzz)) ddpsidzz = psi_spline%deriv(R, z, 0, 2) * psia + if (present(ddpsidrz)) ddpsidrz = psi_spline%deriv(R, z, 1, 1) * psia else - if(present(psinv)) psinv=-1.0_wp_ - if(present(dpsidr)) dpsidr=0.0_wp_ - if(present(dpsidz)) dpsidz=0.0_wp_ - if(present(ddpsidrr)) ddpsidrr=0.0_wp_ - if(present(ddpsidzz)) ddpsidzz=0.0_wp_ - if(present(ddpsidrz)) ddpsidrz=0.0_wp_ + if (present(psi)) psi = -1 + if (present(dpsidr)) dpsidr = 0 + if (present(dpsidz)) dpsidz = 0 + if (present(ddpsidrr)) ddpsidrr = 0 + if (present(ddpsidzz)) ddpsidzz = 0 + if (present(ddpsidrz)) ddpsidrz = 0 end if end subroutine equinum_psi + subroutine equinum_fpol(psi, fpol, dfpol) + ! Computes the poloidal current function F(ψ) + ! (and its derivative) using the numerical equilibrium - subroutine sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc) - use dierckx, only : fpbisp implicit none -! arguments - real(wp_), intent(in) :: rpsim,zpsim - integer, intent(in) :: nur,nuz - real(wp_), intent(out) :: derpsi - real(wp_), dimension(:) :: cc -! local variables - integer, dimension(1) :: iwrkr,iwrkz - real(wp_), dimension(1) :: rrs,zzs,ffspl - real(wp_), dimension(1,ksplp) :: wrkr - real(wp_), dimension(1,ksplp) :: wrkz - rrs(1)=rpsim - zzs(1)=zpsim - call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc,kspl-nur,kspl-nuz, & - rrs,1,zzs,1,ffspl,wrkr,wrkz,iwrkr,iwrkz) - derpsi=ffspl(1)*psia - end subroutine sub_derpsi + ! subroutine arguments + real(wp_), intent(in) :: psi + real(wp_), intent(out) :: fpol + real(wp_), intent(out), optional :: dfpol - - - subroutine equinum_fpol(psinv,fpolv,dfpv) - use dierckx, only : splev,splder - implicit none -! arguments - real(wp_), intent(in) :: psinv - real(wp_), intent(out) :: fpolv - real(wp_), intent(out), optional :: dfpv -! local variables - integer :: ier - real(wp_), dimension(1) :: rrs,ffspl - real(wp_), dimension(nsf) :: wrkfd -! - if(psinv.le.1.0_wp_.and.psinv.ge.0.0_wp_) then - rrs(1)=psinv - call splev(tfp,nsf,cfp,3,rrs,ffspl,1,ier) - fpolv=ffspl(1) - if(present(dfpv)) then - call splder(tfp,nsf,cfp,3,1,rrs,ffspl,1,wrkfd,ier) - dfpv=ffspl(1)/psia - end if + if(psi <= 1 .and. psi >= 0) then + fpol = fpol_spline%eval(psi) + if (present(dfpol)) dfpol = fpol_spline%deriv(psi) / psia else - fpolv=fpolas - if (present(dfpv)) dfpv=0._wp_ + fpol = fpolas + if (present(dfpol)) dfpol = 0 end if end subroutine equinum_fpol - - subroutine equian(rrm,zzm,psinv,fpolv,dfpv,dpsidr,dpsidz, & - ddpsidrr,ddpsidzz,ddpsidrz) - use const_and_precisions, only : wp_ + subroutine equian(R, z, psi, fpolv, dfpv, dpsidr, dpsidz, & + ddpsidrr, ddpsidzz, ddpsidrz) + ! Computes all MHD equilibrium parameters from a simple analytical model implicit none -! arguments - real(wp_), intent(in) :: rrm,zzm - real(wp_), intent(out), optional :: psinv,fpolv,dfpv,dpsidr,dpsidz, & - ddpsidrr,ddpsidzz,ddpsidrz -! local variables - real(wp_) :: cst,dpsidrp,d2psidrp,dqq,qq,rn,rpm,snt,rhop,rhot,btaxqq -! simple model for equilibrium: large aspect ratio -! outside plasma: analytical continuation, not solution Maxwell eqs + ! subroutine arguments + real(wp_), intent(in) :: R,z + real(wp_), intent(out), optional :: psi, fpolv, dfpv, dpsidr, dpsidz, & + ddpsidrr, ddpsidzz, ddpsidrz - rpm=sqrt((rrm-rmaxis)**2+(zzm-zmaxis)**2) !!! rpm==rho_tor[m], rn=rho_tor_norm - rn=rpm/aminor + ! local variables + real(wp_) :: cst, dpsidrp, d2psidrp, dqq, qq, & + rn, rpm, snt, rhop, rhot, btaxqq - snt=0.0_wp_ - cst=1.0_wp_ - if (rpm > 0.0_wp_) then - snt=(zzm-zmaxis)/rpm - cst=(rrm-rmaxis)/rpm + ! simple model for equilibrium: large aspect ratio + ! outside plasma: analytical continuation, not solution Maxwell eqs + + ! Note: rpm is ρ_t in meters, rn is ρ_t normalised + rpm = hypot(R - rmaxis, z - zmaxis) + rn = rpm/aminor + + snt = 0 + cst = 1 + if (rpm > 0) then + snt = (z - zmaxis)/rpm + cst = (R - rmaxis)/rpm end if - if (present(psinv)) then + if (present(psi)) then rhot=rn - if(rn <= 1.0_wp_) then - rhop=frhopol(rhot) - psinv=rhop**2 + if(rn <= 1) then + rhop = frhopol(rhot) + psi = rhop**2 else - psinv=1.0_wp_+btaxis/(2.0_wp_*psia*qa)*(rpm**2-aminor**2) - rhop=sqrt(psinv) + psi = 1 + btaxis / (2*psia*model%q1) * (rpm**2 - aminor**2) + rhop = sqrt(psi) end if end if - if(rn <= 1.0_wp_) then - qq=q0+(qa-q0)*rn**alq - btaxqq=btaxis/qq - dpsidrp=btaxqq*rpm - dqq=alq*(qa-q0)*rn**(alq-1.0_wp_) - d2psidrp=btaxqq*(1.0_wp_-rn*dqq/qq) + if(rn <= 1) then + qq = model%q0 + (model%q1 - model%q0) * rn**model%lq + btaxqq = btaxis/qq + dpsidrp = btaxqq*rpm + dqq = model%lq * (model%q1 - model%q0) * rn**(model%lq - 1) + d2psidrp=btaxqq*(1 - rn*dqq/qq) else - btaxqq=btaxis/qa - dpsidrp=btaxqq*rpm - d2psidrp=btaxqq + btaxqq = btaxis / model%q1 + dpsidrp = btaxqq * rpm + d2psidrp = btaxqq end if - if(present(fpolv)) fpolv=btaxis*rmaxis - if(present(dfpv)) dfpv=0.0_wp_ + if(present(fpolv)) fpolv = btaxis * rmaxis + if(present(dfpv)) dfpv = 0 - if(present(dpsidr)) dpsidr=dpsidrp*cst - if(present(dpsidz)) dpsidz=dpsidrp*snt - if(present(ddpsidrr)) ddpsidrr=btaxqq*snt**2+cst**2*d2psidrp - if(present(ddpsidrz)) ddpsidrz=cst*snt*(d2psidrp-btaxqq) - if(present(ddpsidzz)) ddpsidzz=btaxqq*cst**2+snt**2*d2psidrp + if(present(dpsidr)) dpsidr = dpsidrp*cst + if(present(dpsidz)) dpsidz = dpsidrp*snt + if(present(ddpsidrr)) ddpsidrr = btaxqq*snt**2 + cst**2*d2psidrp + if(present(ddpsidrz)) ddpsidrz = cst * snt * (d2psidrp - btaxqq) + if(present(ddpsidzz)) ddpsidzz = btaxqq * cst**2 + snt**2 * d2psidrp end subroutine equian - function frhopol(rhot) - use utils, only : locate - use simplespline, only : spli - implicit none -! arguments - real(wp_), intent(in) :: rhot - real(wp_) :: frhopol -! local variables - integer :: i - real(wp_) :: dr + ! Converts from toroidal (ρ_t) to poloidal (ρ_p) normalised radius - call locate(rhotr,nrho,rhot,i) - i=min(max(1,i),nrho-1) - dr=rhot-rhotr(i) - frhopol=spli(crhop,nrho,i,dr) + implicit none + + ! function arguments + real(wp_), intent(in) :: rhot + real(wp_) :: frhopol + + frhopol = rhop_spline%eval(rhot) end function frhopol - function frhopolv(rhot) + ! Vector variant of `frhopol` use utils, only : locate - use simplespline, only : spli + implicit none -! arguments - real(wp_), dimension(:), intent(in) :: rhot - real(wp_), dimension(size(rhot)) :: frhopolv -! local variables - integer :: i,i0,j + + ! function arguments + real(wp_), intent(in) :: rhot(:) + real(wp_) :: frhopolv(size(rhot)) + + ! local variables + integer :: i, i0, j real(wp_) :: dr - i0=1 - do j=1,size(rhot) - call locate(rhotr(i0:),nrho-i0+1,rhot(j),i) - i=min(max(1,i),nrho-i0)+i0-1 - dr=rhot(j)-rhotr(i) - frhopolv(j)=spli(crhop,nrho,i,dr) - i0=i + i0 = 1 + do j = 1, size(rhot) + call locate(rhop_spline%data(i0:), rhop_spline%ndata-i0+1, rhot(j), i) + i = min(max(1,i), rhop_spline%ndata - i0) + i0 - 1 + dr = rhot(j) - rhop_spline%data(i) + frhopolv(j) = rhop_spline%raw_eval(i, dr) + i0 = i end do end function frhopolv - function frhotor(rhop) - use utils, only : locate - use simplespline, only : spli - implicit none -! arguments - real(wp_), intent(in) :: rhop - real(wp_) :: frhotor -! local variables - integer :: i - real(wp_) :: dr + ! Converts from poloidal (ρ_p) to toroidal (ρ_t) normalised radius - call locate(rhopr,nrho,rhop,i) - i=min(max(1,i),nrho-1) - dr=rhop-rhopr(i) - frhotor=spli(crhot,nrho,i,dr) + implicit none + + ! function arguments + real(wp_), intent(in) :: rhop + real(wp_) :: frhotor + + frhotor = rhot_spline%eval(rhop) end function frhotor - function fq(psin) - use const_and_precisions, only : wp_ + ! Computes the safety factor q as a function of the + ! normalised poloidal flux ψ use gray_params, only : iequil - use simplespline, only :spli - use utils, only : locate + implicit none -! arguments + + ! function arguments real(wp_), intent(in) :: psin real(wp_) :: fq -! local variables - integer :: i - real(wp_) :: dps,rn - if (iequil<2) then - rn=frhotor(sqrt(psin)) - fq=q0+(qa-q0)*rn**alq + ! local variables + real(wp_) :: rn + + if (iequil < 2) then + ! q = q₀ + (q₁ - q₀)ρ^l + rn = frhotor(sqrt(psin)) + fq = model%q0 + (model%q1 - model%q0) * rn**model%lq else - call locate(psinr,nq,psin,i) - i=min(max(1,i),nq-1) - dps=psin-psinr(i) - fq=spli(cq,nq,i,dps) + fq = q_spline%eval(psin) end if end function fq - - subroutine bfield(rpsim,zpsim,bphi,br,bz) + subroutine bfield(rpsim, zpsim, bphi, br, bz) + ! Computes the magnetic field as a function of + ! (R, z) in cylindrical coordinates use gray_params, only : iequil + implicit none -! arguments - real(wp_), intent(in) :: rpsim,zpsim + + ! subroutine arguments + real(wp_), intent(in) :: rpsim, zpsim real(wp_), intent(out), optional :: bphi,br,bz -! local variables + + ! local variables real(wp_) :: psin,fpol if (iequil < 2) then call equian(rpsim,zpsim,fpolv=bphi,dpsidr=bz,dpsidz=br) if (present(bphi)) bphi=bphi/rpsim else - call equinum_psi(rpsim,zpsim,psinv=bphi,dpsidr=bz,dpsidz=br) + call equinum_psi(rpsim,zpsim,psi=bphi,dpsidr=bz,dpsidz=br) if (present(bphi)) then psin=bphi call equinum_fpol(psin,fpol) @@ -1083,34 +1031,52 @@ contains end subroutine bfield + function tor_curr(rpsim,zpsim) result(jphi) + ! Computes the toroidal current Jφ as a function of (R, z) + use const_and_precisions, only : ccj=>mu0inv + use gray_params, only : iequil - subroutine tor_curr(rpsim,zpsim,ajphi) - use const_and_precisions, only : wp_,ccj=>mu0inv - use gray_params, only : iequil implicit none -! arguments - real(wp_) :: rpsim,zpsim,ajphi -! local variables + + ! function arguments + real(wp_), intent(in) :: rpsim, zpsim + real(wp_) :: jphi + + ! local variables real(wp_) :: bzz,dbvcdc13,dbvcdc31 real(wp_) :: dpsidr,ddpsidrr,ddpsidzz if(iequil < 2) then - call equian(rpsim,zpsim,dpsidr=dpsidr, & + call equian(rpsim,zpsim,dpsidr=dpsidr, & ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz) else - call equinum_psi(rpsim,zpsim,dpsidr=dpsidr, & + call equinum_psi(rpsim,zpsim,dpsidr=dpsidr, & ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz) end if bzz= dpsidr/rpsim dbvcdc13=-ddpsidzz/rpsim dbvcdc31= ddpsidrr/rpsim-bzz/rpsim - ajphi=ccj*(dbvcdc13-dbvcdc31) - end subroutine tor_curr + jphi=ccj*(dbvcdc13-dbvcdc31) + end function tor_curr + function tor_curr_psi(psin) result(jphi) + ! Computes the toroidal current Jφ as a function of ψ + + implicit none + + ! function arguments + real(wp_), intent(in) :: psin + real(wp_) :: jphi + + ! local variables + real(wp_) :: r1, r2 + call psi_raxis(psin, r1, r2) + jphi = tor_curr(r2, zmaxis) + end function tor_curr_psi + subroutine psi_raxis(psin,r1,r2) - use const_and_precisions, only : wp_ use gray_params, only : iequil use dierckx, only : profil, sproota use logger, only : log_error @@ -1127,7 +1093,7 @@ contains integer :: iopt,ier,m real(wp_) :: zc,val real(wp_), dimension(mest) :: zeroc - real(wp_), dimension(nsr) :: czc + real(wp_), dimension(psi_spline%nknots_x) :: czc character(64) :: msg if (iequil < 2) then @@ -1137,35 +1103,26 @@ contains else iopt=1 zc=zmaxis - call profil(iopt,tr,nsr,tz,nsz,cceq,kspl,kspl,zc,nsr,czc,ier) + call profil(iopt, psi_spline%knots_x, psi_spline%nknots_x, & + psi_spline%knots_y, psi_spline%nknots_y, & + psi_spline%coeffs, kspl, kspl, zc, & + psi_spline%nknots_x, czc, ier) if (ier > 0) then write (msg, '("profil failed with error ",g0)') ier call log_error(msg, mod='equilibrium', proc='psi_raxis') end if val=psin*psiant+psinop - call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) + call sproota(val, psi_spline%knots_x, psi_spline%nknots_x, & + czc, zeroc, mest, m, ier) r1=zeroc(1) r2=zeroc(2) end if end subroutine psi_raxis - - subroutine tor_curr_psi(psin,ajphi) - use const_and_precisions, only : wp_ - implicit none -! arguments - real(wp_) :: psin,ajphi -! local variables - real(wp_) :: r1,r2 - call psi_raxis(psin,r1,r2) - call tor_curr(r2,zmaxis,ajphi) - end subroutine tor_curr_psi - - - subroutine points_ox(rz,zz,rf,zf,psinvf,info) + ! Finds the location of the O,X points use const_and_precisions, only : comp_eps use minpack, only : hybrj1 use logger, only : log_error, log_debug @@ -1203,7 +1160,6 @@ contains end subroutine points_ox - subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag) use logger, only : log_error @@ -1238,7 +1194,6 @@ contains end subroutine fcnox - subroutine points_tgo(rz,zz,rf,zf,psin0,info) use const_and_precisions, only : comp_eps use minpack, only : hybrj1mv @@ -1278,10 +1233,8 @@ contains end - subroutine fcntgo(n,x,f0,fvec,fjac,ldfjac,iflag) - use const_and_precisions, only : wp_ - use logger, only : log_error + use logger, only : log_error implicit none @@ -1314,45 +1267,15 @@ contains end subroutine fcntgo - subroutine unset_equil_spline + ! Unsets the splines global variables, see the top of this file. implicit none - if (allocated(tr)) deallocate(tr) - if (allocated(tz)) deallocate(tz) - if (allocated(tfp)) deallocate(tfp) - if (allocated(cfp)) deallocate(cfp) - if (allocated(cceq)) deallocate(cceq) - if (allocated(cceq01)) deallocate(cceq01) - if (allocated(cceq10)) deallocate(cceq10) - if (allocated(cceq02)) deallocate(cceq02) - if (allocated(cceq20)) deallocate(cceq20) - if (allocated(cceq11)) deallocate(cceq11) - nsr = 0 - nsz = 0 - nsf = 0 + call fpol_spline%deinit + call psi_spline%deinit + call q_spline%deinit + call rhop_spline%deinit + call rhot_spline%deinit end subroutine unset_equil_spline - - - subroutine unset_q - implicit none - - if (allocated(psinr)) deallocate(psinr) - if (allocated(cq)) deallocate(cq) - nq = 0 - end subroutine unset_q - - - - subroutine unset_rho_spline - implicit none - - if(allocated(rhopr)) deallocate(rhopr) - if(allocated(rhotr)) deallocate(rhotr) - if(allocated(crhop)) deallocate(crhop) - if(allocated(crhot)) deallocate(crhot) - nrho=0 - end subroutine unset_rho_spline - end module equilibrium diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 556157e..746cfa1 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -2187,9 +2187,9 @@ bb: do subroutine print_prof ! Prints the (input) plasma kinetic profiles (unit 98) - use equilibrium, only : psinr, nq, fq, frhotor, tor_curr_psi - use coreprofiles, only : density, temp - use units, only : uprfin, unit_active + use equilibrium, only : q_spline, fq, frhotor, tor_curr_psi + use coreprofiles, only : density, temp + use units, only : uprfin, unit_active implicit none @@ -2198,19 +2198,19 @@ bb: do ! local variables integer :: i - real(wp_) :: psin, rhot, ajphi, dens, ddens + real(wp_) :: psin, rhot, jphi, dens, ddens if (.not. unit_active(uprfin)) return write (uprfin, *) '#psi rhot ne Te q Jphi' - do i=1,nq - psin = psinr(i) + do i = 1, q_spline%ndata + psin = q_spline%data(i) rhot = frhotor(sqrt(psin)) call density(psin, dens, ddens) - call tor_curr_psi(max(eps, psin), ajphi) + jphi = tor_curr_psi(max(eps, psin)) write (uprfin, '(12(1x,e12.5))') & - psin, rhot, dens, temp(psin), fq(psin), ajphi*1.e-6_wp_ + psin, rhot, dens, temp(psin), fq(psin), jphi*1.e-6_wp_ end do end subroutine print_prof @@ -2218,8 +2218,8 @@ bb: do subroutine print_bres(bres) ! Prints the EC resonance surface table (unit 70) - use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield, nq - use units, only : ubres, unit_active + use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield, q_spline + use units, only : ubres, unit_active implicit none @@ -2234,14 +2234,14 @@ bb: do integer, dimension(10) :: ncpts real(wp_) :: dr,dz,btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb real(wp_), dimension(icmx) :: rrcb,zzcb - real(wp_) :: rv(nq), zv(nq) - real(wp_), dimension(nq,nq) :: btotal + real(wp_) :: rv(q_spline%ndata), zv(q_spline%ndata) + real(wp_), dimension(q_spline%ndata,q_spline%ndata) :: btotal if (.not. unit_active(ubres)) return - dr = (rmxm - rmnm)/(nq - 1) - dz = (zmxm - zmnm)/(nq - 1) - do j=1,nq + dr = (rmxm - rmnm)/(q_spline%ndata - 1) + dz = (zmxm - zmnm)/(q_spline%ndata - 1) + do j=1,q_spline%ndata rv(j) = rmnm + dr*(j - 1) zv(j) = zmnm + dz*(j - 1) end do @@ -2249,9 +2249,9 @@ bb: do ! Btotal on psi grid btmx = -1.0e30_wp_ btmn = 1.0e30_wp_ - do k=1,nq + do k = 1, q_spline%ndata zzk = zv(k) - do j=1,nq + do j = 1, q_spline%ndata rrj = rv(j) call bfield(rrj, zzk, bbphi, bbr, bbz) btotal(j,k) = sqrt(bbr**2 + bbz**2 + bbphi**2) @@ -2262,14 +2262,14 @@ bb: do ! compute Btot=Bres/n with n=1,5 write (ubres, *) '#i Btot R z' - do n=1,5 + do n = 1, 5 bbb = bres/dble(n) if (bbb >= btmn .and. bbb <= btmx) then nconts = size(ncpts) nctot = size(rrcb) - call cniteq(rv, zv, btotal, nq, nq, bbb, & + call cniteq(rv, zv, btotal, q_spline%ndata, q_spline%ndata, bbb, & nconts, ncpts, nctot, rrcb, zzcb) - do j=1,nctot + do j = 1, nctot write (ubres, '(i6,12(1x,e12.5))') j, bbb, rrcb(j), zzcb(j) end do end if @@ -2284,7 +2284,7 @@ bb: do use gray_params, only : iequil use equilibrium, only : rmnm, rmxm, zmnm, zmxm, equian, equinum_psi, & - equinum_fpol, nq + equinum_fpol, q_spline use coreprofiles, only : density, temp use units, only : umaps, unit_active @@ -2297,28 +2297,28 @@ bb: do integer :: j,k real(wp_) :: dr, dz, zk, rj, bphi, br, bz, btot, & psin, ne, dne, te, xg, yg, anpl - real(wp_), dimension(nq) :: r, z + real(wp_), dimension(q_spline%ndata) :: r, z if (.not. unit_active(umaps)) return - dr = (rmxm-rmnm)/(nq - 1) - dz = (zmxm-zmnm)/(nq - 1) - do j=1,nq + dr = (rmxm-rmnm)/(q_spline%ndata - 1) + dz = (zmxm-zmnm)/(q_spline%ndata - 1) + do j=1,q_spline%ndata r(j) = rmnm + dr*(j - 1) z(j) = zmnm + dz*(j - 1) end do write (umaps, *) '#R z psin Br Bphi Bz Btot ne Te X Y Npl' - do j=1,nq + do j = 1, q_spline%ndata rj = r(j) anpl = anpl0 * r0/rj - do k=1,nq + do k = 1, q_spline%ndata zk = z(k) if (iequil < 2) then - call equian(rj, zk, psinv=psin, fpolv=bphi, dpsidr=bz, dpsidz=br) + call equian(rj, zk, psi=psin, fpolv=bphi, dpsidr=bz, dpsidz=br) else - call equinum_psi(rj, zk, psinv=psin, dpsidr=bz, dpsidz=br) - call equinum_fpol(psin, fpolv=bphi) + call equinum_psi(rj, zk, psi=psin, dpsidr=bz, dpsidz=br) + call equinum_fpol(psin, fpol=bphi) end if br = -br/rj bphi = bphi/rj @@ -2340,7 +2340,7 @@ bb: do subroutine print_surfq(qval) ! Print constant ψ surfaces for a given `q` value - use equilibrium, only : psinr, nq, fq, frhotor, & + use equilibrium, only : q_spline, fq, frhotor, & rmaxis, zmaxis, zbsup, zbinf use magsurf_data, only : contours_psi, npoints, print_contour use utils, only : locate, intlin @@ -2355,23 +2355,23 @@ bb: do integer :: i1,i real(wp_) :: rup,zup,rlw,zlw,rhot,psival real(wp_), dimension(npoints) :: rcn,zcn - real(wp_), dimension(nq) :: qpsi + real(wp_), dimension(q_spline%ndata) :: qpsi character(256) :: msg ! for log messages formatting - ! build q profile on psin grid - do i=1,nq - qpsi(i) = fq(psinr(i)) + ! build the q profile on the ψ grid + do i = 1, q_spline%ndata + qpsi(i) = fq(q_spline%data(i)) end do ! locate ψ surface for q=qval call log_info('constant ψ surfaces for:', & mod='gray_core', proc='print_surfq') - do i=1,size(qval) + do i=1, size(qval) ! FIXME: check for non monotonous q profile - call locate(abs(qpsi),nq,qval(i),i1) - if (i1>0 .and. i1 0 .and. i1 < q_spline%ndata) then + call intlin(abs(qpsi(i1)), q_spline%data(i1), abs(qpsi(i1+1)), & + q_spline%data(i1+1), qval(i), psival) rup = rmaxis rlw = rmaxis zup = (zbsup + zmaxis)/2.0_wp_ diff --git a/src/gray_jetto1beam.f90 b/src/gray_jetto1beam.f90 index 76bc3f3..d080dbe 100644 --- a/src/gray_jetto1beam.f90 +++ b/src/gray_jetto1beam.f90 @@ -127,13 +127,11 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, & ! Free memory free_memory: block - use equilibrium, only : unset_equil_spline, unset_rho_spline, unset_q + use equilibrium, only : unset_equil_spline use coreprofiles, only : unset_profiles_spline ! Unset global variables of the `equilibrium` module call unset_equil_spline - call unset_rho_spline - call unset_q ! Unset global variables of the `coreprofiles` module call unset_profiles_spline diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90 index 6c2662e..5509d83 100644 --- a/src/magsurf_data.f90 +++ b/src/magsurf_data.f90 @@ -1,5 +1,7 @@ module magsurf_data use const_and_precisions, only : wp_ + use splines, only : spline_simple + implicit none integer, save :: npsi, npoints !# sup mag, # punti per sup @@ -14,12 +16,11 @@ module magsurf_data real(wp_), dimension(:), allocatable, save :: tjp,tlm,ch,ch01 real(wp_), dimension(:,:), allocatable, save :: rcon,zcon - real(wp_), dimension(:,:), allocatable, save :: cdadrhot,cdvdrhot - real(wp_), dimension(:,:), allocatable, save :: cvol,crri,crbav,cbmx,cbmn,carea,cfc - real(wp_), dimension(:,:), allocatable, save :: crhotq - real(wp_), dimension(:,:), allocatable, save :: cratja,cratjb,cratjpl - + type(spline_simple), save :: cvol, crri, crbav, cbmx, cbmn, carea, cfc + type(spline_simple), save :: crhotq + type(spline_simple), save :: cratja, cratjb, cratjpl + type(spline_simple), save :: cdadrhot, cdvdrhot contains @@ -58,46 +59,43 @@ contains allocate(pstab(npsi), & rhot_eq(npsi),rhotqv(npsi),bav(npsi),bmxpsi(npsi),bmnpsi(npsi),varea(npsi), & vvol(npsi),vcurrp(npsi),vajphiav(npsi),qqv(npsi),ffc(npsi),vratja(npsi), & - vratjb(npsi),rpstab(npsi),rri(npsi),rbav(npsi),cdadrhot(npsi,4), & - cdvdrhot(npsi,4),cbmx(npsi,4),cbmn(npsi,4),crbav(npsi,4),cvol(npsi,4), & - crri(npsi,4),carea(npsi,4),cfc(npsi,4),crhotq(npsi,4),cratjpl(npsi,4), & - cratja(npsi,4),cratjb(npsi,4)) - end subroutine alloc_surfvec + vratjb(npsi),rpstab(npsi),rri(npsi),rbav(npsi)) + end subroutine alloc_surfvec subroutine dealloc_surfvec implicit none call dealloc_cnt - if(allocated(pstab)) deallocate(pstab) + if(allocated(pstab)) deallocate(pstab) if(allocated(rhot_eq)) deallocate(rhot_eq) - if(allocated(rhotqv)) deallocate(rhotqv) - if(allocated(bav)) deallocate(bav) - if(allocated(bmxpsi)) deallocate(bmxpsi) - if(allocated(bmnpsi)) deallocate(bmnpsi) - if(allocated(varea)) deallocate(varea) - if(allocated(vvol)) deallocate(vvol) - if(allocated(vcurrp)) deallocate(vcurrp) - if(allocated(vajphiav)) deallocate(vajphiav) - if(allocated(qqv)) deallocate(qqv) - if(allocated(ffc)) deallocate(ffc) - if(allocated(vratja)) deallocate(vratja) - if(allocated(vratjb)) deallocate(vratjb) - if(allocated(rpstab)) deallocate(rpstab) - if(allocated(rri)) deallocate(rri) - if(allocated(rbav)) deallocate(rbav) - if(allocated(cdadrhot)) deallocate(cdadrhot) - if(allocated(cdvdrhot)) deallocate(cdvdrhot) - if(allocated(cbmx)) deallocate(cbmx) - if(allocated(cbmn)) deallocate(cbmn) - if(allocated(crbav)) deallocate(crbav) - if(allocated(cvol)) deallocate(cvol) - if(allocated(crri)) deallocate(crri) - if(allocated(carea)) deallocate(carea) - if(allocated(cfc)) deallocate(cfc) - if(allocated(crhotq)) deallocate(crhotq) - if(allocated(cratjpl)) deallocate(cratjpl) - if(allocated(cratja)) deallocate(cratja) - if(allocated(cratjb)) deallocate(cratjb) - if(allocated(tjp)) deallocate(tjp,tlm,ch) + if(allocated(rhotqv)) deallocate(rhotqv) + if(allocated(bav)) deallocate(bav) + if(allocated(bmxpsi)) deallocate(bmxpsi) + if(allocated(bmnpsi)) deallocate(bmnpsi) + if(allocated(varea)) deallocate(varea) + if(allocated(vvol)) deallocate(vvol) + if(allocated(vcurrp)) deallocate(vcurrp) + if(allocated(vajphiav)) deallocate(vajphiav) + if(allocated(qqv)) deallocate(qqv) + if(allocated(ffc)) deallocate(ffc) + if(allocated(vratja)) deallocate(vratja) + if(allocated(vratjb)) deallocate(vratjb) + if(allocated(rpstab)) deallocate(rpstab) + if(allocated(rri)) deallocate(rri) + if(allocated(rbav)) deallocate(rbav) + if(allocated(tjp)) deallocate(tjp,tlm,ch) + + call cvol%deinit + call crbav%deinit + call crri%deinit + call cbmx%deinit + call cbmn%deinit + call cratja%deinit + call cratjb%deinit + call cratjpl%deinit + call carea%deinit + call cfc%deinit + call cdadrhot%deinit + call cdvdrhot%deinit end subroutine dealloc_surfvec @@ -106,8 +104,7 @@ contains use gray_params, only : iequil use equilibrium, only : btrcen,btaxis,rmaxis,zmaxis,phitedge,zbsup,zbinf, & equian,equinum_psi,bfield,frhotor,fq,tor_curr - use simplespline, only : difcs - use dierckx, only : regrid,coeff_parder + use dierckx, only : regrid,coeff_parder implicit none ! local constants @@ -232,7 +229,7 @@ contains bmmx=-1.0e+30_wp_ bmmn=1.0e+30_wp_ - call tor_curr(rctemp(1),zctemp(1),ajphi0) + ajphi0 = tor_curr(rctemp(1),zctemp(1)) call bfield(rctemp(1),zctemp(1),bphi,br=brr,bz=bzz) fpolv=bphi*rctemp(1) btot0=sqrt(bphi**2+brr**2+bzz**2) @@ -260,7 +257,7 @@ contains rpsim=rctemp(inc1) zpsim=zctemp(inc1) call bfield(rpsim,zpsim,br=brr,bz=bzz) - call tor_curr(rpsim,zpsim,ajphi) + ajphi = tor_curr(rpsim,zpsim) bphi=fpolv/rpsim btot=sqrt(bphi**2+brr**2+bzz**2) bpoloid=sqrt(brr**2+bzz**2) @@ -372,36 +369,22 @@ contains rpstab(npsi)=1.0_wp_ pstab(npsi)=1.0_wp_ -! spline coefficients of area,vol,rbav,rri,bmxpsi,bmnpsi,fc,dadrhot,dvdrhot,ratioJs -! used for computations of dP/dV and J_cd - iopt=0 - call difcs(rpstab,vvol,npsi,iopt,cvol,ier) - iopt=0 - call difcs(rpstab,rbav,npsi,iopt,crbav,ier) - iopt=0 - call difcs(rpstab,rri,npsi,iopt,crri,ier) - iopt=0 - call difcs(rpstab,bmxpsi,npsi,iopt,cbmx,ier) - iopt=0 - call difcs(rpstab,bmnpsi,npsi,iopt,cbmn,ier) - iopt=0 - call difcs(rpstab,vratja,npsi,iopt,cratja,ier) - iopt=0 - call difcs(rpstab,vratjb,npsi,iopt,cratjb,ier) - iopt=0 - call difcs(rpstab,vratjpl,npsi,iopt,cratjpl,ier) - iopt=0 - call difcs(rpstab,varea,npsi,iopt,carea,ier) - iopt=0 - call difcs(rpstab,ffc,npsi,iopt,cfc,ier) - iopt=0 - call difcs(rpstab,dadrhotv,npsi,iopt,cdadrhot,ier) - iopt=0 - call difcs(rpstab,dvdrhotv,npsi,iopt,cdvdrhot,ier) -! iopt=0 -! call difcs(rpstab,qqv,npsi,iopt,cqq,ier) + ! spline coefficients of area,vol,rbav,rri,bmxpsi,bmnpsi,fc,dadrhot,dvdrhot,ratioJs + ! used for computations of dP/dV and J_cd + call cvol%init(rpstab, vvol, npsi) + call crbav%init(rpstab, rbav, npsi) + call crri%init(rpstab, rri, npsi) + call cbmx%init(rpstab, bmxpsi, npsi) + call cbmn%init(rpstab, bmnpsi, npsi) + call cratja%init(rpstab, vratja, npsi) + call cratjb%init(rpstab, vratjb, npsi) + call cratjpl%init(rpstab, vratjpl, npsi) + call carea%init(rpstab, varea, npsi) + call cfc%init(rpstab, ffc, npsi) + call cdadrhot%init(rpstab, dadrhotv, npsi) + call cdvdrhot%init(rpstab, dvdrhotv, npsi) -! spline interpolation of H(lambda,rhop) and dH/dlambda + ! spline interpolation of H(lambda,rhop) and dH/dlambda iopt=0 s=0.0_wp_ call regrid(iopt,npsi,rpstab,nlam,alam,ffhlam,zero,one,zero,one, & @@ -430,37 +413,30 @@ contains subroutine fluxval(rhop,area,vol,dervol,dadrhot,dvdrhot, & rri,rbav,bmn,bmx,fc,ratja,ratjb,ratjpl) use const_and_precisions, only : wp_ - use utils, only : locate - use simplespline, only :spli,splid implicit none -! arguments + + ! subroutine arguments real(wp_), intent(in) :: rhop - real(wp_), intent(out), optional :: vol,area,rri,rbav,dervol,bmn,bmx,fc, & - ratja,ratjb,ratjpl,dadrhot,dvdrhot -! local variables - integer :: ip - real(wp_) :: drh + real(wp_), intent(out), optional :: & + vol, area, rri, rbav, dervol, bmn, bmx, fc, & + ratja, ratjb, ratjpl, dadrhot, dvdrhot - call locate(rpstab,npsi,rhop,ip) - ip=min(max(1,ip),npsi-1) - drh=rhop-rpstab(ip) + if (present(area)) area = carea%eval(rhop) + if (present(vol)) vol = cvol%eval(rhop) - if (present(area)) area=spli(carea,npsi,ip,drh) - if (present(vol)) vol=spli(cvol,npsi,ip,drh) + if (present(dervol)) dervol = cvol%deriv(rhop) + if (present(dadrhot)) dadrhot = cdadrhot%eval(rhop) + if (present(dvdrhot)) dvdrhot = cdvdrhot%eval(rhop) - if (present(dervol)) dervol=splid(cvol,npsi,ip,drh) - if (present(dadrhot)) dadrhot=spli(cdadrhot,npsi,ip,drh) - if (present(dvdrhot)) dvdrhot=spli(cdvdrhot,npsi,ip,drh) + if (present(rri)) rri = crri%eval(rhop) + if (present(rbav)) rbav = crbav%eval(rhop) + if (present(bmn)) bmn = cbmn%eval(rhop) + if (present(bmx)) bmx = cbmx%eval(rhop) + if (present(fc)) fc = cfc%eval(rhop) - if (present(rri)) rri=spli(crri,npsi,ip,drh) - if (present(rbav)) rbav=spli(crbav,npsi,ip,drh) - if (present(bmn)) bmn=spli(cbmn,npsi,ip,drh) - if (present(bmx)) bmx=spli(cbmx,npsi,ip,drh) - if (present(fc)) fc=spli(cfc,npsi,ip,drh) - - if (present(ratja)) ratja=spli(cratja,npsi,ip,drh) - if (present(ratjb)) ratjb=spli(cratjb,npsi,ip,drh) - if (present(ratjpl)) ratjpl=spli(cratjpl,npsi,ip,drh) + if (present(ratja)) ratja = cratja%eval(rhop) + if (present(ratjb)) ratjb = cratjb%eval(rhop) + if (present(ratjpl)) ratjpl = cratjpl%eval(rhop) end subroutine fluxval @@ -470,22 +446,26 @@ contains use const_and_precisions, only : wp_,pi use gray_params, only : iequil use logger, only : log_warning - use dierckx, only : profil,sproota - use equilibrium, only : rmaxis,zmaxis,aminor,frhotor,tr,nsr,tz,nsz,cceq, & - kspl,psiant,psinop,points_tgo - use limiter, only : rwallm + use dierckx, only : profil, sproota + use equilibrium, only : rmaxis, zmaxis, aminor, frhotor, psi_spline, & + kspl, psiant, psinop, points_tgo + use limiter, only : rwallm + implicit none -! local constants + + ! local constants integer, parameter :: mest=4 -! arguments + + ! subroutine arguments real(wp_), intent(in) :: h real(wp_), dimension(:), intent(out) :: rcn,zcn real(wp_), intent(inout) :: rup,zup,rlw,zlw -! local variables + + ! local variables integer :: npoints,np,info,ic,ier,iopt,m real(wp_) :: ra,rb,za,zb,rn,th,zc,val real(wp_), dimension(mest) :: zeroc - real(wp_), dimension(nsr) :: czc + real(wp_), dimension(psi_spline%nknots_x) :: czc npoints=size(rcn) np=(npoints-1)/2 @@ -517,7 +497,10 @@ contains do ic=2,np zc=zlw+(zup-zlw)*(1.0_wp_-cos(th*(ic-1)))/2.0_wp_ iopt=1 - call profil(iopt,tr,nsr,tz,nsz,cceq,kspl,kspl,zc,nsr,czc,ier) + call profil(iopt, psi_spline%knots_x, psi_spline%nknots_x, & + psi_spline%knots_y, psi_spline%nknots_y, & + psi_spline%coeffs, kspl, kspl, zc, & + psi_spline%nknots_x, czc, ier) if (ier > 0) then block character(256) :: msg @@ -527,7 +510,8 @@ contains end block end if val=h*psiant+psinop - call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) + call sproota(val, psi_spline%knots_x, psi_spline%nknots_x, & + czc, zeroc, mest, m, ier) if (zeroc(1).gt.rwallm) then rcn(ic)=zeroc(1) rcn(npoints+1-ic)=zeroc(2) diff --git a/src/main.f90 b/src/main.f90 index bd00a63..804c270 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -193,7 +193,7 @@ contains ! or an analytical description) and initialises the respective ! GRAY parameters and data. use equilibrium, only : read_equil_an, read_eqdsk, change_cocos, & - set_equil_an, set_equil_spline, eq_scal + set_equil_an, set_equil_spline, scale_equil use logger, only : log_debug implicit none @@ -223,7 +223,7 @@ contains end if ! Rescale B, I and/or force their signs - call eq_scal(params%equilibrium, data%equilibrium) + call scale_equil(params%equilibrium, data%equilibrium) ! Set global variables (for splines) if (params%equilibrium%iequil < 2) then @@ -239,7 +239,7 @@ contains subroutine deinit_equilibrium(data) ! Free all memory allocated by the init_equilibrium subroutine. use gray_params, only : equilibrium_data - use equilibrium, only : unset_equil_spline, unset_rho_spline, unset_q + use equilibrium, only : unset_equil_spline implicit none @@ -254,8 +254,6 @@ contains ! Unset global variables of the `equilibrium` module call unset_equil_spline - call unset_rho_spline - call unset_q end subroutine deinit_equilibrium diff --git a/src/simplespline.f90 b/src/simplespline.f90 deleted file mode 100644 index b1c5b68..0000000 --- a/src/simplespline.f90 +++ /dev/null @@ -1,273 +0,0 @@ -module simplespline - - use const_and_precisions, only : wp_ - implicit none - -contains - -function spli(cspli,n,k,dx) - implicit none - integer, intent(in) :: n, k - real(wp_), intent(in) :: cspli(n,4), dx - real(wp_) :: spli - spli=cspli(k,1)+dx*(cspli(k,2)+dx*(cspli(k,3)+dx*cspli(k,4))) -end function spli - -function splid(cspli,n,k,dx) - implicit none - integer, intent(in) :: n, k - real(wp_), intent(in) :: cspli(n,4), dx - real(wp_) :: splid - splid=cspli(k,2)+dx*(2.0_wp_*cspli(k,3)+3.0_wp_*dx*cspli(k,4)) -end function splid - -subroutine difcs(x,y,n,iopt,c,ier) - implicit none - integer, intent(in) :: n, iopt - real(wp_), intent(in) :: x(n), y(n) - real(wp_), intent(inout) :: c(n*4) - integer :: ier - integer :: jmp,iol,ioh,i,ii,j,j1,j2,j3 - real(wp_) :: xb,xc,ya,yb,h,a,r,dya,dyb,dy2 - jmp =1 - if (n <= 1) return -! -! initialization -! - xc =x(1) - yb =y(1) - h =0.0_wp_ - a =0.0_wp_ - r =0.0_wp_ - dyb =0.0_wp_ -! -! iol=0 - given derivative at first point -! ioh=0 - given derivative at last point -! - iol =iopt-1 - ioh =iopt-2 - if (ioh == 1) then - iol =0 - ioh =0 - end if - dy2 =c(2) -! -! form the system of linear equations -! and eliminate subsequentially -! - j =1 - do i=1,n - j2 =n+i - j3 =j2+n - a =h*(2.0_wp_-a) - dya =dyb+h*r - if (i>=n) then -! -! set derivative dy2 at last point -! - dyb =dy2 - h =0.0_wp_ - if (ioh/=0) then - dyb =dya - goto 13 - end if - else - j =j+jmp - xb =xc - xc =x(j) - h =xc-xb -! -! ii=0 - increasing abscissae -! ii=1 - decreasing abscissae -! - ii =0 - if (h==0) return - if (h<0) ii =1 - ya =yb - yb =y(j) - dyb =(yb-ya)/h - if (i<=1) then - j1 =ii - if (iol/=0) goto 13 - dya =c(1) - end if - end if - if (j1-ii /= 0) return - a =1.0_wp_/(h+h+a) - 13 continue - r =a*(dyb-dya) - c(j3)=r - a =h*a - c(j2)=a - c(i) =dyb - end do -! -! back substitution of the system of linear equations -! and computation of the other coefficients -! - a =1.0_wp_ - j1 =j3+n+ii-ii*n - i =n - do iol=1,n - xb =x(j) - h =xc-xb - xc =xb - a =a+h - yb =r - r =c(j3)-r*c(j2) - ya =r+r - c(j3)=ya+r - c(j2)=c(i)-h*(ya+yb) - c(j1)=(yb-r)/a - c(i) =y(j) - a =0.0_wp_ - j =j-jmp - i =i-1 - j2 =j2-1 - j3 =j3-1 - j1 =j3+n+ii - end do - ier =0 -end subroutine difcs - -subroutine difcsn(xx,yy,nmx,n,iopt,cc,ier) -! -! same as difcs but with dimension(xx,yy) = nmx > n -! - implicit none - integer, intent(in) :: nmx, n, iopt - real(wp_), intent(in) :: xx(nmx), yy(nmx) - real(wp_), intent(inout) :: cc(nmx,4) - integer :: ier - integer :: jmp,iol,ioh,i,ii,j,j1,j2,j3 - real(wp_) :: x(n),y(n),c(n*4),xb,xc,ya,yb,h,a,r,dya,dyb,dy2 -! - do i=1,n - x(i)=xx(i) - y(i)=yy(i) - end do - ii=0 - do j=1,4 - do i=1,n - ii=ii+1 - c(ii)=cc(i,j) - end do - end do -! - jmp =1 - if (n>1) then -! -! initialization -! - xc =x(1) - yb =y(1) - h =0.0_wp_ - a =0.0_wp_ - r =0.0_wp_ - dyb =0.0_wp_ -! -! iol=0 - given derivative at first point -! ioh=0 - given derivative at last point -! - iol =iopt-1 - ioh =iopt-2 - if (ioh==1) then - iol =0 - ioh =0 - end if - dy2 =c(2) -! -! form the system of linear equations -! and eliminate subsequentially -! - j =1 - do i=1,n - j2 =n+i - j3 =j2+n - a =h*(2.0_wp_-a) - dya =dyb+h*r - if (i>=n) then -! -! set derivative dy2 at last point -! - dyb =dy2 - h =0.0_wp_ - if (ioh/=0) then - dyb =dya - goto 13 - end if - else - j =j+jmp - xb =xc - xc =x(j) - h =xc-xb -! -! ii=0 - increasing abscissae -! ii=1 - decreasing abscissae -! - ii =0 - if (h==0) goto 16 - if (h<0) ii =1 - ya =yb - yb =y(j) - dyb =(yb-ya)/h - if (i<=1) then - j1 =ii - if (iol/=0) goto 13 - dya =c(1) - end if - end if - if (j1/=ii) goto 16 - a =1.0_wp_/(h+h+a) - 13 continue - r =a*(dyb-dya) - c(j3)=r - a =h*a - c(j2)=a - c(i) =dyb - end do -! -! back substitution of the system of linear equations -! and computation of the other coefficients -! - a =1.0_wp_ - j1 =j3+n+ii-ii*n - i =n - do iol=1,n - xb =x(j) - h =xc-xb - xc =xb - a =a+h - yb =r - r =c(j3)-r*c(j2) - ya =r+r - c(j3)=ya+r - c(j2)=c(i)-h*(ya+yb) - c(j1)=(yb-r)/a - c(i) =y(j) - a =0.0_wp_ - j =j-jmp - i =i-1 - j2 =j2-1 - j3 =j3-1 - j1 =j3+n+ii - end do - ier =0 - end if -! - 16 continue - ii=0 - do j=1,4 - do i=1,nmx - if(i<=n) then - ii=ii+1 - cc(i,j)=c(ii) - else - cc(i,j)=0.0_wp_ - end if - end do - end do -! -end subroutine difcsn - -end module simplespline \ No newline at end of file diff --git a/src/splines.f90 b/src/splines.f90 new file mode 100644 index 0000000..7174323 --- /dev/null +++ b/src/splines.f90 @@ -0,0 +1,546 @@ +! This module provides a high-level interface for creating and evaluating +! several kind of splines: +! +! `spline_simple` is a simple interpolating cubic spline, +! `spline_1d` and `spline_2d` are wrapper around the DIERCKX cubic splines. +module splines + use const_and_precisions, only : wp_ + + implicit none + + ! A 1D interpolating cubic spline + type spline_simple + integer :: ndata ! Number of data points + real(wp_), allocatable :: data(:) ! Data points (ndata) + real(wp_), allocatable :: coeffs(:,:) ! Spline coefficients (ndata, 4) + contains + procedure :: init => spline_simple_init + procedure :: deinit => spline_simple_deinit + procedure :: eval => spline_simple_eval + procedure :: raw_eval => spline_simple_raw_eval + procedure :: deriv => spline_simple_deriv + end type + + ! A 1D smoothing/interpolating cubic spline + type spline_1d + integer :: nknots ! Number of spline knots + real(wp_), allocatable :: knots(:) ! Knots positions + real(wp_), allocatable :: coeffs(:) ! B-spline coefficients + contains + procedure :: init => spline_1d_init + procedure :: deinit => spline_1d_deinit + procedure :: eval => spline_1d_eval + procedure :: deriv => spline_1d_deriv + end type + + ! A 2D smoothing/interpolating cubic spline s(x, y) + type spline_2d + integer :: nknots_x ! Number of x knots + integer :: nknots_y ! Number of y knots + real(wp_), allocatable :: knots_x(:) ! Knots x positions + real(wp_), allocatable :: knots_y(:) ! Knots y positions + real(wp_), allocatable :: coeffs(:) ! B-spline coefficients + + ! B-spline coefficients of the partial derivatives + type(pointer), allocatable :: partial(:,:) + + contains + procedure :: init => spline_2d_init + procedure :: deinit => spline_2d_deinit + procedure :: eval => spline_2d_eval + procedure :: init_deriv => spline_2d_init_deriv + procedure :: deriv => spline_2d_deriv + end type + + ! Wrapper to store pointers in an array + type pointer + real(wp_), pointer :: ptr(:) => null() + end type + + private + public spline_simple, spline_1d, spline_2d + +contains + + subroutine spline_simple_init(self, x, y, n) + ! Initialises the spline + + implicit none + + ! subroutine arguments + class(spline_simple), intent(inout) :: self + integer, intent(in) :: n + real(wp_), dimension(n), intent(in) :: x, y + + call self%deinit + self%ndata = n + allocate(self%data(n)) + allocate(self%coeffs(n, 4)) + + self%data = x + call spline_simple_coeffs(x, y, n, self%coeffs) + end subroutine spline_simple_init + + + subroutine spline_simple_deinit(self) + ! Deinitialises a simple_spline + implicit none + + ! subroutine arguments + class(spline_simple), intent(inout) :: self + + if (allocated(self%data)) deallocate(self%data) + if (allocated(self%coeffs)) deallocate(self%coeffs) + self%ndata = 0 + end subroutine spline_simple_deinit + + + function spline_simple_eval(self, x) result(y) + ! Evaluates the spline at x + use utils, only : locate + + implicit none + + ! subroutine arguments + class(spline_simple), intent(in) :: self + real(wp_), intent(in) :: x + real(wp_) :: y + + ! local variables + integer :: i + real(wp_) :: dx + + call locate(self%data, self%ndata, x, i) + i = min(max(1, i), self%ndata - 1) + dx = x - self%data(i) + y = self%raw_eval(i, dx) + end function spline_simple_eval + + + function spline_simple_raw_eval(self, i, dx) result(y) + ! Evaluates the i-th polynomial of the spline at dx + + implicit none + + ! subroutine arguments + class(spline_simple), intent(in) :: self + integer, intent(in) :: i + real(wp_), intent(in) :: dx + real(wp_) :: y + + y = self%coeffs(i,1) + dx*(self%coeffs(i,2) & + + dx*(self%coeffs(i,3) + dx*self%coeffs(i,4))) + end function spline_simple_raw_eval + + + function spline_simple_deriv(self, x) result(y) + ! Computes the derivative of the spline at x + use utils, only : locate + + implicit none + + ! subroutine arguments + class(spline_simple), intent(in) :: self + real(wp_), intent(in) :: x + real(wp_) :: y + + ! local variables + integer :: i + real(wp_) :: dx + + call locate(self%data, self%ndata, x, i) + i = min(max(1, i), self%ndata - 1) + dx = x - self%data(i) + y = self%coeffs(i,2) + dx*(2*self%coeffs(i,3) + 3*dx*self%coeffs(i,4)) + end function spline_simple_deriv + + + subroutine spline_simple_coeffs(x, y, n, c) + ! Computes the cubic coefficients of all n polynomials + implicit none + + ! subroutine arguments + integer, intent(in) :: n + real(wp_), intent(in) :: x(n), y(n) + real(wp_), intent(inout) :: c(n*4) + + ! local variables + integer :: jmp, i, ii, j, j1, j2, j3, k + real(wp_) :: xb, xc, ya, yb, h, a, r, dya, dyb, dy2 + jmp = 1 + + if (n <= 1) return + + ! initialisation + xc = x(1) + yb = y(1) + h = 0 + a = 0 + r = 0 + dyb = 0 + dy2 = c(2) + + ! form the system of linear equations + ! and eliminate subsequentially + j = 1 + do i = 1, n + j2 = n + i + j3 = j2 + n + a = h*(2 - a) + dya = dyb + h*r + if (i >= n) then + ! set derivative dy2 at last point + dyb = dy2 + h = 0 + dyb = dya + goto 13 + else + j = j+jmp + xb = xc + xc = x(j) + h = xc-xb + + ! ii=0 - increasing abscissae + ! ii=1 - decreasing abscissae + ii = 0 + if (h == 0) return + if (h < 0) ii = 1 + ya = yb + yb = y(j) + dyb = (yb - ya)/h + + if (i <= 1) then + j1 = ii + goto 13 + end if + end if + + if (j1-ii /= 0) return + a = 1 / (2*h + a) + 13 continue + r = a*(dyb - dya) + c(j3) = r + a = h*a + c(j2) = a + c(i) = dyb + end do + + ! back substitution of the system of linear equations + ! and computation of the other coefficients + a = 1 + j1 = j3+n+ii-ii*n + i = n + do k = 1, n + xb = x(j) + h = xc - xb + xc = xb + a = a+h + yb = r + r = c(j3)-r*c(j2) + ya = 2*r + c(j3) = ya + r + c(j2) = c(i) - h*(ya+yb) + c(j1) = (yb - r)/a + c(i) = y(j) + a = 0 + j = j-jmp + i = i-1 + j2 = j2-1 + j3 = j3-1 + j1 = j3+n+ii + end do + end subroutine spline_simple_coeffs + + + subroutine spline_1d_init(self, x, y, n, range, weights, tension, err) + ! Initialises a spline_1d. + ! Takes: + ! x: x data points + ! y: y data points + ! n: number of data points + ! range: interpolation range as [x_min, x_max] + ! weights: factors weighting the data points (default: all 1) + ! tension: parameter controlling the amount of smoothing (default: 0) + ! Returns: + ! err: error code of `curfit` + use dierckx, only : curfit + + implicit none + + ! subroutine arguments + class(spline_1d), intent(inout) :: self + real(wp_), intent(in) :: x(n) + real(wp_), intent(in) :: y(n) + integer, intent(in) :: n + real(wp_), intent(in) :: range(2) + real(wp_), intent(in), optional :: weights(n) + real(wp_), intent(in), optional :: tension + integer, intent(out), optional :: err + + ! local variables + integer :: nknots_est ! over-estimate of the number of knots + real(wp_) :: residuals ! sum of the residuals + integer :: work_int(n + 4) ! integer working space + real(wp_) :: work_real(4*n + 16*(n + 4)) ! real working space + + ! default values + integer :: err_def + real(wp_) :: weights_def(n), tension_def + + weights_def = 1 + tension_def = 0 + if (present(weights)) weights_def = weights + if (present(tension)) tension_def = tension + + ! clear memory, if necessary + call self%deinit + + ! allocate the spline arrays + nknots_est = n + 4 + allocate(self%knots(nknots_est), self%coeffs(nknots_est)) + + call curfit(0, n, x, y, weights_def, range(1), range(2), 3, tension_def, & + nknots_est, self%nknots, self%knots, self%coeffs, residuals, & + work_real, size(work_real), work_int, err_def) + + if (present(err)) err = err_def + end subroutine spline_1d_init + + + subroutine spline_1d_deinit(self) + ! Deinitialises a spline_1d + implicit none + + class(spline_1d), intent(inout) :: self + + if (allocated(self%knots)) deallocate(self%knots) + if (allocated(self%coeffs)) deallocate(self%coeffs) + self%nknots = 0 + end subroutine spline_1d_deinit + + + function spline_1d_eval(self, x) result(y) + ! Evaluates the spline at x + use dierckx, only : splev + + implicit none + + ! subroutine arguments + class(spline_1d), intent(in) :: self + real(wp_), intent(in) :: x + real(wp_) :: y + + ! local variables + integer :: err + real(wp_) :: yv(1) ! because splev returns a vector + + call splev(self%knots, self%nknots, self%coeffs, 3, [x], yv, 1, err) + y = yv(1) + end function spline_1d_eval + + + function spline_1d_deriv(self, x, order) result(y) + ! Evaluates the spline n-th order derivative at x + use dierckx, only : splder + + implicit none + + ! subroutine arguments + class(spline_1d), intent(in) :: self + real(wp_), intent(in) :: x + integer, intent(in), optional :: order + real(wp_) :: y + + ! local variables + integer :: err, n + real(wp_) :: yv(1) ! because splev returns a vector + real(wp_) :: work(self%nknots) ! working space array + + n = 1 + if (present(order)) n = order + + call splder(self%knots, self%nknots, self%coeffs, & + 3, n, [x], yv, 1, work, err) + y = yv(1) + end function spline_1d_deriv + + + subroutine spline_2d_init(self, x, y, z, nx, ny, range, tension, err) + ! Initialises a spline_2d. + ! Takes: + ! x, y: data points on a regular grid + ! z: data points of z(x, y) + ! n: number of data points + ! range: interpolation range as [x_min, x_max, y_min, y_max] + ! weights: factors weighting the data points (default: all 1) + ! tension: parameter controlling the amount of smoothing (default: 0) + ! Returns: + ! err: error code of `curfit` + use dierckx, only : regrid + + implicit none + + ! subroutine arguments + class(spline_2d), intent(inout) :: self + real(wp_), intent(in) :: x(nx) + real(wp_), intent(in) :: y(ny) + real(wp_), intent(in) :: z(nx * ny) + integer, intent(in) :: nx, ny + real(wp_), intent(in) :: range(4) + real(wp_), intent(in), optional :: tension + integer, intent(out), optional :: err + + ! local variables + integer :: nknots_x_est ! over-estimate of the number of knots + integer :: nknots_y_est ! + real(wp_) :: residuals ! sum of the residuals + + ! working space arrays + integer :: work_int(2*(ny + nx) + 11) + real(wp_) :: work_real(15*(nx + ny) + ny*(nx + 4) + 92 + max(ny, nx + 4)) + + ! default values + integer :: err_def + real(wp_) :: tension_def + + tension_def = 0 + if (present(tension)) tension_def = tension + + ! clear memory, if necessary + call self%deinit + + ! allocate the spline arrays + nknots_x_est = nx + 4 + nknots_y_est = ny + 4 + allocate(self%knots_x(nknots_x_est), self%knots_y(nknots_y_est)) + allocate(self%coeffs(nx * ny)) + + call regrid(0, nx, x, ny, y, z, range(1), range(2), range(3), range(4), & + 3, 3, tension_def, nknots_x_est, nknots_y_est, & + self%nknots_x, self%knots_x, self%nknots_y, self%knots_y, & + self%coeffs, residuals, work_real, size(work_real), & + work_int, size(work_int), err_def) + + if (present(err)) err = err_def + end subroutine spline_2d_init + + + subroutine spline_2d_deinit(self) + ! Deinitialises a spline_2d + implicit none + + class(spline_2d), intent(inout) :: self + + if (allocated(self%knots_x)) deallocate(self%knots_x) + if (allocated(self%knots_y)) deallocate(self%knots_y) + if (allocated(self%coeffs)) deallocate(self%coeffs) + + ! Note: partial derivatives coeff. are pointers + if (allocated(self%partial)) then + deallocate(self%partial(1, 0)%ptr) + deallocate(self%partial(0, 1)%ptr) + deallocate(self%partial(1, 1)%ptr) + deallocate(self%partial(2, 0)%ptr) + deallocate(self%partial(0, 2)%ptr) + deallocate(self%partial) + end if + + self%nknots_x = 0 + self%nknots_y = 0 + end subroutine spline_2d_deinit + + + function spline_2d_eval(self, x, y) result(z) + ! Evaluates the spline at (x, y) + use dierckx, only : fpbisp + + implicit none + + ! subroutine arguments + class(spline_2d), intent(in) :: self + real(wp_), intent(in) :: x, y + real(wp_) :: z + + ! local variables + real(wp_) :: zv(1) ! because fpbisp returns a vector + + integer :: work_int(8) ! integer working space + real(wp_) :: work_real(8) ! real working space + + ! Note: see https://netlib.org/dierckx/bispev.f for + ! this apparently nonsensical invocation + call fpbisp(self%knots_x, self%nknots_x, & + self%knots_y, self%nknots_y, & + self%coeffs, 3, 3, [x], 1, [y], 1, & + zv, work_real(1), work_real(5), & + work_int(1), work_int(2)) + z = zv(1) + end function spline_2d_eval + + + subroutine spline_2d_init_deriv(self, p, q, n, m) + ! Computes the spline coefficients of n-th partial derivative + ! w.r.t x and m-th partial derivative w.r.t y on a grid of + ! p×q points. + ! + ! Note: for simplicity, only up to second-order is supported. + use dierckx, only : coeff_parder + + implicit none + + ! subroutine arguments + class(spline_2d), intent(inout) :: self + integer, intent(in) :: p, q ! grid dimensions + integer, intent(in) :: n, m ! derivative order + + ! local variables + integer :: coeff_size + integer :: err + + ! coeff. array (actually, the working space) size + coeff_size = p*(4 - n) + q*(4 - m) + p*q + + ! allocate slots for storing the derivatives (first call only) + if (.not. allocated(self%partial)) allocate(self%partial(0:2, 0:2)) + + ! allocate the coefficients array + allocate(self%partial(n, m)%ptr(coeff_size)) + + ! compute the coefficients + call coeff_parder(self%knots_x, self%nknots_x, & + self%knots_y, self%nknots_y, & + self%coeffs, 3, 3, n, m, & + self%partial(n, m)%ptr, coeff_size, err) + end subroutine spline_2d_init_deriv + + + function spline_2d_deriv(self, x, y, n, m) result(z) + ! Evaluates the spline n-th partial derivative w.r.t x + ! and m-th partial derivative w.r.t y at (x, y) + ! + ! Note: the coefficients of the derivative must have been + ! initialised with init_deriv before calling this method. + use dierckx, only : fpbisp + + implicit none + + ! subroutine arguments + class(spline_2d), intent(in) :: self + real(wp_), intent(in) :: x, y + integer, intent(in) :: n, m + real(wp_) :: z + + ! local variables + real(wp_), dimension(1) :: zv ! because splev returns a vector + integer, dimension(1) :: work_int_x, work_int_y ! integer working space + real(wp_), dimension(1,4) :: work_real_x, work_real_y ! real working space + + call fpbisp(self%knots_x(1 + n), self%nknots_x - 2*n, & + self%knots_y(1 + m), self%nknots_y - 2*m, & + self%partial(n, m)%ptr, & + 3 - n, 3 - m, [x], 1, [y], 1, zv, & + work_real_x, work_real_y, work_int_x, work_int_y) + z = zv(1) + end function spline_2d_deriv + +end module splines