src/splines.f90: add spline_2d%init_nonreg

This adds a proper subroutine to initialise a spline_2d given
non-regular data using surfit.

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_) :: psinoptmp, psinxptmp
-    real(wp_) :: rbmin, rbmax, rbinf, rbsup, r1, z1
+    real(wp_) :: rbinf, rbsup, r1, z1
     real(wp_) :: psiant, psinop
-    real(wp_),  dimension(size(data%psinr)) :: rhotn
+    real(wp_) :: rhotn(size(data%psinr))
-    real(wp_),  dimension(:),  allocatable :: rv1d, zv1d, fvpsi, wf
+    real(wp_), allocatable :: rv1d(:), zv1d(:), fvpsi(:)
-    integer :: ixploc, info, i, j, ij
+    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 psi_spline%init_nonreg(rv1d, zv1d, fvpsi, tension=params%ssplps, &
-        call scatterspl(rv1d, zv1d, fvpsi, wf, nrz, kspl, tension, &
+                                    range=[rmnm, rmxm, zmnm, zmxm], err=err)
-                        rmnm, rmxm, zmnm, zmxm,                    &
-                        psi_spline%knots_x, psi_spline%nknots_x,   &
-                        psi_spline%knots_y, psi_spline%nknots_y,   &
-                        psi_spline%coeffs, 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, &
+          call psi_spline%init_nonreg(rv1d, zv1d, fvpsi, tension=zero, &
-                          rmnm, rmxm, zmnm, zmxm,                    &
+                                      range=[rmnm, rmxm, zmnm, zmxm], err=err)
-                          psi_spline%knots_x, psi_spline%nknots_x,   &
-                          psi_spline%knots_y, psi_spline%nknots_y,   &
-                          psi_spline%coeffs, 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))
+    block
-    wf(1:npsi-1) = 1
+      real(wp_) :: w(npsi)
-    wf(npsi)     = 1.0e2_wp_
+      w(1:npsi-1) = 1
+      w(npsi)     = 1.0e2_wp_
       call fpol_spline%init(data%psinr, data%fpol, npsi, range=[zero, one], &
-                          weights=wf, tension=params%ssplf)
+                            weights=w, tension=params%ssplf)
-    deallocate(wf)
+    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
+    call find_ox_point(data%rax,data%zax,rmaxis,zmaxis,psinoptmp)
-    zax0=data%zax
-    call find_ox_point(rax0,zax0,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 (params%ixp)
-
-    select case (ixploc)
       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

src/splines.f90

@@ -52,6 +52,7 @@ module splines
 
     contains
       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
@@ -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