diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index f26bee7..646cc09 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -408,12 +408,12 @@ contains ! local variables 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_) :: psinoptmp, psinxptmp + real(wp_) :: rbinf, rbsup, r1, z1 real(wp_) :: psiant, psinop - real(wp_), dimension(size(data%psinr)) :: rhotn - real(wp_), dimension(:), allocatable :: rv1d, zv1d, fvpsi, wf - integer :: ixploc, info, i, j, ij + real(wp_) :: rhotn(size(data%psinr)) + real(wp_), allocatable :: rv1d(:), zv1d(:), fvpsi(:) + integer :: i, j, ij character(256) :: msg ! for log messages formatting ! compute array sizes @@ -461,7 +461,7 @@ contains end do ! store valid data - allocate(rv1d(nrz),zv1d(nrz),fvpsi(nrz),wf(nrz)) + allocate(rv1d(nrz),zv1d(nrz),fvpsi(nrz)) ij=0 do j=1,nz do i=1,nr @@ -474,7 +474,6 @@ contains rv1d(ij)=data%rv(i) zv1d(ij)=data%zv(j) fvpsi(ij)=data%psin(i,j) - wf(ij)=1.0d0 end do end do @@ -482,26 +481,18 @@ contains ! use reduced number of knots to limit memory comsumption ? 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, err) + call psi_spline%init_nonreg(rv1d, zv1d, fvpsi, tension=params%ssplps, & + range=[rmnm, rmxm, zmnm, zmxm], err=err) ! if failed, re-fit with an interpolating spline (zero tension) if(err == -1) then err = 0 - 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, err) + call psi_spline%init_nonreg(rv1d, zv1d, fvpsi, tension=zero, & + range=[rmnm, rmxm, zmnm, zmxm], err=err) end if - deallocate(rv1d, zv1d, wf, fvpsi) + deallocate(rv1d, zv1d, fvpsi) ! reset nrz to the total number of grid points for next allocations nrz = nr*nz @@ -549,12 +540,13 @@ contains ! Spline interpolation of F(ψ) ! give a small weight to the last point - allocate(wf(npsi)) - 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) + block + real(wp_) :: w(npsi) + w(1:npsi-1) = 1 + w(npsi) = 1.0e2_wp_ + call fpol_spline%init(data%psinr, data%fpol, npsi, range=[zero, one], & + weights=w, tension=params%ssplf) + end block ! set vacuum value used outside 0 ≤ ψ ≤ 1 range fpolas = fpol_spline%eval(data%psinr(npsi)) @@ -575,29 +567,22 @@ contains call vmaxmini(data%boundary%z, nbnd, zbinf, zbsup, ibinf, ibsup) rbinf = data%boundary%R(ibinf) rbsup = data%boundary%R(ibsup) - call vmaxmin(data%boundary%R, nbnd, rbmin, rbmax) else zbinf=data%zv(2) zbsup=data%zv(nz-1) rbinf=data%rv((nr+1)/2) rbsup=rbinf - rbmin=data%rv(2) - rbmax=data%rv(nr-1) end if ! Search for exact location of the magnetic axis - rax0=data%rax - zax0=data%zax - call find_ox_point(rax0,zax0,rmaxis,zmaxis,psinoptmp) + call find_ox_point(data%rax,data%zax,rmaxis,zmaxis,psinoptmp) write (msg, '("O-point found:", 3(x,a,"=",g0.3))') & 'r', rmaxis, 'z', zmaxis, 'ψ', psinoptmp call log_info(msg, mod='equilibrium', proc='set_equil_spline') ! Search for X-point - ixploc = params%ixp - - select case (ixploc) + select case (params%ixp) case (X_AT_BOTTOM) call find_ox_point(rbinf,zbinf,r1,z1,psinxptmp) if(psinxptmp/=-1.0_wp_) then @@ -610,8 +595,6 @@ contains psiant=psinxptmp-psinop call find_htg_point(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one) zbsup=z1 - else - ixploc=0 end if case (X_AT_TOP) @@ -626,8 +609,6 @@ contains psiant=psinxptmp-psinop call find_htg_point(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one) zbinf=z1 - else - ixploc=0 end if case (X_IS_MISSING) @@ -725,57 +706,6 @@ contains end subroutine rescale_boundary - subroutine scatterspl(x,y,z,w,n,kspl,sspl,xmin,xmax,ymin,ymax, & - tx,nknt_x,ty,nknt_y,coeff,ierr) - ! 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 - - ! subroutine arguments - integer, intent(in) :: n - real(wp_), dimension(n), intent(in) :: x, y, z - real(wp_), dimension(n), intent(in) :: w - integer, intent(in) :: kspl - real(wp_), intent(in) :: sspl - real(wp_), intent(in) :: xmin, xmax, ymin, ymax - real(wp_), dimension(nknt_x), intent(inout) :: tx - real(wp_), dimension(nknt_y), intent(inout) :: ty - integer, intent(inout) :: nknt_x, nknt_y - real(wp_), dimension(nknt_x*nknt_y), intent(out) :: coeff - integer, intent(out) :: ierr - - ! local variables - integer :: iopt - real(wp_) :: resid - integer :: u,v,km,ne,b1,b2,lwrk1,lwrk2,kwrk,nxest,nyest - real(wp_), dimension(:), allocatable :: wrk1, wrk2 - integer, dimension(:), allocatable :: iwrk - - nxest=nknt_x - nyest=nknt_y - ne = max(nxest,nyest) - - km = kspl+1 - u = nxest-km - v = nyest-km - b1 = kspl*min(u,v)+kspl+1 - b2 = (kspl+1)*min(u,v)+1 - lwrk1 = u*v*(2+b1+b2)+2*(u+v+km*(n+ne)+ne-2*kspl)+b2+1 - lwrk2 = u*v*(b2+1)+b2 - kwrk = n+(nknt_x-2*kspl-1)*(nknt_y-2*kspl-1) - allocate(wrk1(lwrk1),wrk2(lwrk2),iwrk(kwrk)) - - iopt=0 - call surfit(iopt,n,x,y,z,w,xmin,xmax,ymin,ymax,kspl,kspl, & - sspl,nxest,nyest,ne,comp_eps,nknt_x,tx,nknt_y,ty, & - coeff,resid,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ierr) - - deallocate(wrk1,wrk2,iwrk) - - end subroutine scatterspl - - subroutine setqphi_num(psinq,q,psia,rhotn) ! Computes the spline of the safety factor q(ψ) use const_and_precisions, only : pi diff --git a/src/splines.f90 b/src/splines.f90 index 94f3825..e1175d5 100644 --- a/src/splines.f90 +++ b/src/splines.f90 @@ -51,12 +51,13 @@ module splines 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 - procedure :: transform => spline_2d_transform + procedure :: init => spline_2d_init + procedure :: init_nonreg => spline_2d_init_nonreg + procedure :: deinit => spline_2d_deinit + procedure :: eval => spline_2d_eval + procedure :: init_deriv => spline_2d_init_deriv + procedure :: deriv => spline_2d_deriv + procedure :: transform => spline_2d_transform end type ! A simple 1D linear interpolation @@ -370,10 +371,9 @@ contains ! 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` + ! err: error code of `regrid` use dierckx, only : regrid ! subroutine arguments @@ -421,6 +421,68 @@ contains end subroutine spline_2d_init + subroutine spline_2d_init_nonreg(self, x, y, z, weights, tension, range, err) + ! Initialises a spline_2d for non-regular data. + ! Takes: + ! x, y: data points (not necessarily 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 `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 + + ! subroutine arguments + class(spline_2d), intent(inout) :: self + real(wp_), intent(in) :: x(:), y(:), z(:) + real(wp_), intent(in) :: range(4) + real(wp_), intent(in), optional :: weights(:) + real(wp_), intent(in), optional :: tension + integer, intent(out), optional :: err + + ! local variables + real(wp_) :: residuals + real(wp_) :: tension_def, weights_def(size(z)) + integer :: err_def + + tension_def = 0 + if (present(tension)) tension_def = tension + + weights_def = 1 + if (present(weights)) weights_def = weights + + associate(n => size(z), & + u => self%nknots_x - 4, & + v => self%nknots_y - 4, & + ne => max(self%nknots_x, self%nknots_y), & + b1 => 3*min(self%nknots_x, self%nknots_y), & + b2 => 4*min(self%nknots_x, self%nknots_y) - 3) + block + ! working-space arrays + real(wp_) :: wrk1(u*v*(2+b1+b2) + 2*(u+v+4*(n+ne)+ne) + b2 - 11) + real(wp_) :: wrk2(u*v*(1+b2) + b2) + integer :: iwrk(n + (self%nknots_x-7)*(self%nknots_y-7)) + + call surfit(0, size(z), x, y, z, weights_def, & + range(1), range(2), range(3), range(4), & + 3, 3, tension_def, & + self%nknots_x, self%nknots_x, ne, comp_eps, & + self%nknots_x, self%knots_x, & + self%nknots_y, self%knots_y, & + self%coeffs, residuals, wrk1, size(wrk1), & + wrk2, size(wrk2), iwrk, size(iwrk), err_def) + end block + end associate + + if (present(err)) err = err_def + end subroutine spline_2d_init_nonreg + + subroutine spline_2d_deinit(self) ! Deinitialises a spline_2d