module interpolation_m use screen_io_m, only : get_stderr use precision_m, only: FP use descriptors_m, only : ERR_UNHANDLED use error_handling_m, only: handle_error, error_info_t use status_codes_m, only : PARALLAX_ERR_INTERPOL implicit none public interpol_coeffs public interpol1d private interpol1d_linear private interpol1d_cubic contains subroutine interpol_coeffs(intorder,x,y,coeffs) !****************************************************************************** ! computes coefficients of interpolation such that quantity at position x,y can be computed as ! ! u(x,y) = sum_{i=1,intorder+1,j=1,intorder+1} coeffs(i,j) u(i,j) ! ! i,j is Cartesian ginterpolation grid with x_i=0...intorder, y_j=0...intorder ! x,y is coordinate of to be interpolated point with respect to southwest point of interpolation stamp maesured in grid distance ! ! INPUT: ! intorder: order of interpolation ! x: x-coordinate w.r.t southwest node point ! y: y-coordinate w.r.t southwest node point ! ! OUTPUT: ! coeffs: coefficients of interpolating polynomial ! !****************************************************************************** implicit none integer,intent(in)::intorder real(kind=FP),intent(in)::x,y real(kind=FP),intent(out),dimension(intorder+1,intorder+1)::coeffs integer::nads,i,j,k,info,lwork integer,dimension(intorder+1)::ipiv real(kind=FP),dimension(intorder+1,intorder+1)::V real(kind=FP),dimension(intorder+1)::work real(kind=FP),dimension(intorder+1)::g,h if (intorder.le.0) then write(get_stderr(), *) & 'error(interpol_coeffs): intorder must be larger 0 but is', intorder ERROR STOP ERR_UNHANDLED endif nads=intorder+1 ! Build Vandermonde Matrix do j=1,nads do i=1,nads V(i,j)=1.0_FP*(i-1)**(j-1) enddo enddo ! invert Vandermonde Matrix ! critical statements are necessary when multithreaded ! lapack/blas libraries are used as called from parallel region !$omp critical #ifdef DOUBLE_PREC call dgetrf(nads,nads,V,nads,ipiv,info) #else call sgetrf(nads,nads,V,nads,ipiv,info) #endif !$omp end critical if (info.ne.0) then write(get_stderr(), *) 'error(interpol_coeffs): dgetrf, info=', info ERROR STOP ERR_UNHANDLED endif lwork=nads !$omp critical #ifdef DOUBLE_PREC call dgetri(nads,V,nads,ipiv,work,lwork,info) #else call sgetri(nads,V,nads,ipiv,work,lwork,info) #endif !$omp end critical if (info.ne.0) then write(get_stderr(), *) 'error(interpol_coeffs): dgetri, info=', info ERROR STOP ERR_UNHANDLED endif ! interpolation along x (same for all rows) g=0.0_FP do i=1,nads do k=1,nads g(i)=g(i)+x**(k-1)*V(k,i) enddo enddo ! interpolation along y h=0.0_FP do j=1,nads do k=1,nads h(j)=h(j)+y**(k-1)*V(k,j) enddo enddo do i=1,nads do j=1,nads coeffs(i,j)=g(i)*h(j) enddo enddo end subroutine interpol_coeffs real(FP) function interpol1d(intorder, nnodes, dx, ui, x) !! Performs a 1D polynomial interpolation integer, intent(in) :: intorder !! Interpoaltion order, yet only 1 and 3 ar available integer, intent(in) :: nnodes !! Number of nodes of interpolation data real(FP), intent(in) :: dx !! Distance between nodes (equidistant) real(FP), intent(in), dimension(nnodes) :: ui !! Values of interpolation data real(FP), intent(in) :: x if (intorder == 1) then interpol1d = interpol1d_linear(nnodes, dx, ui, x) elseif (intorder == 3) then interpol1d = interpol1d_cubic(nnodes, dx, ui, x) else call handle_error('Only intorder=1, 3 are allowed', & PARALLAX_ERR_INTERPOL, __LINE__, __FILE__, & additional_info=error_info_t('intorder=',[intorder])) endif end function real(FP) function interpol1d_cubic(nnodes, dx, ui, x) !! Performs a 1D cubic interpolation on equally spaced data sets !! Based on Vandermonde matrix with equidistant spacing shifted to !! leftmost points of interpolation stencil. !! V = [1 0 0 0 !! 1 dx 2dx 3dx !! 1 dx^2 (2dx)^2 (3dx^2) !! 1 dx^3 (2dx)^3 (3dx^3)] !! The inverse of this matrix is built explicitly as vndrm_inv integer, intent(in) :: nnodes !! Number of nodes of interpolation data real(FP), intent(in) :: dx !! Distance between nodes (equidistant) real(FP), intent(in), dimension(nnodes) :: ui !! Values of interpolation data real(FP), intent(in) :: x !! Location where interpoaltion data is requested integer, parameter :: intorder = 3 ! Interpolation order integer, parameter :: nint = intorder + 1 integer :: i_left, k real(FP) :: xrel, x_left real(FP), dimension(nint) :: yi, cf real(FP) :: dxinv, dxinv_sqrd, dxinv_cubed real(FP), dimension(nint, nint) :: vndrm_inv real(FP), dimension(nint, nint), parameter :: vndrm_inv_coeffs=reshape(& [1.0_FP, -11.0_FP/6.0_FP, 1.0_FP, -1.0_FP/6.0_FP, & 0.0_FP, 3.0_FP, -5.0_FP/2.0_FP, 1.0_FP/2.0_FP, & 0.0_FP, -3.0_FP/2.0_FP, 2.0_FP, -1.0_FP/2.0_FP, & 0.0_FP, 1.0_FP / 3.0_FP, -1.0_FP/2.0_FP, 1.0_FP/6.0_FP], & [nint, nint]) ! Coefficents in inverse Vandermonde matrix, that are independent of dx if (nnodes < nint) then call handle_error('Number of nodes must be larger than 3', & PARALLAX_ERR_INTERPOL, __LINE__, __FILE__, & additional_info=error_info_t('nnodes=',[nnodes])) endif ! Build inverse of Vandermonde Matrx dxinv = 1.0_FP / dx dxinv_sqrd = dxinv**2 dxinv_cubed = dxinv**3 vndrm_inv(1, 1:4) = vndrm_inv_coeffs(1, 1:4) vndrm_inv(2, 1:4) = vndrm_inv_coeffs(2, 1:4) * dxinv vndrm_inv(3, 1:4) = vndrm_inv_coeffs(3, 1:4) * dxinv_sqrd vndrm_inv(4, 1:4) = vndrm_inv_coeffs(4, 1:4) * dxinv_cubed ! Determine first (left) point of interpolation stamp i_left = floor(x / dx) if (i_left < 1) then i_left = 1 endif if (i_left + intorder > nnodes) then i_left = nnodes - intorder endif x_left = (i_left - 1) * dx ! Determine coefficients of interpolating polynomial do k = 1, nint yi(k) = ui(i_left + k - 1) enddo cf = MatMul(vndrm_inv, yi) ! return solution interpol1d_cubic = 0.0_FP do k = 1, nint xrel = x - x_left interpol1d_cubic = interpol1d_cubic + cf(k) * xrel**(k-1) enddo end function pure real(kind=FP) function interpol1d_linear(nnodes, dx, ui, x) !! Performs 1D linear interpolation/extrapolation !! for set of equally spaced data integer, intent(in) :: nnodes !! Number of nodes real(kind=FP), intent(in) :: dx !! Node spacing real(kind=FP), intent(in), dimension(nnodes) :: ui !! Values at nodes in increasing order real(kind=FP), intent(in) :: x !! Desired interpolation position with respect to first node point integer :: bin ! Bin (an index) into which x falls, ! such that ui(bin) <= interpol1d_linear(x) <= ui(bin+1). bin = floor( x / dx ) + 1 ! If bin < 1 or bin > nnodes - 1, perform extrapolation bin = max( 1, bin ) bin = min( nnodes - 1, bin ) interpol1d_linear = ui(bin) + (x - dx * (bin - 1)) / dx * (ui(bin + 1) - ui(bin)) end function interpol1d_linear end module