rprofile_bspline.F90

!!matvec with matmul !#define __MATVEC_N(y,Mat,Vec) y=MATMUL(Mat,Vec) !#define __MATVEC_T(y,Mat,Vec) y=MATMUL(Vec,Mat) !#define __PMATVEC_N(fy,y,Mat,Vec) y=fyy+MATMUL(Mat,Vec) !#define __PMATVEC_T(fy,y,Mat,Vec) y=fyy+MATMUL(Vec,Mat) !#define __AMATVEC_N(y,fMat,Mat,Vec) y=fMatMATMUL(Mat,Vec) !#define __AMATVEC_T(y,fMat,Mat,Vec) y=fMatMATMUL(Vec,Mat) !#define __PAMATVEC_N(fy,y,fMat,Mat,Vec) y=fyy+fMatMATMUL(Mat,Vec) !#define __PAMATVEC_T(fy,y,fMat,Mat,Vec) y=fyy+fMatMATMUL(Vec,Mat)

!!#define __GENERICMATVEC(NT,fy,y,fMat,Mat,Vec) CALL DGEMV(NT,SIZE(Mat,1),SIZE(Mat,2),fMat,Mat,SIZE(Mat,1),Vec,1,fy,y,1)

!!matmat with matmul !#define __MATMAT_NN(Y,A,B) Y=MATMUL(A,B) !#define __MATMAT_TN(Y,A,B) Y=MATMUL(TRANSPOSE(A),B) !#define __MATMAT_NT(Y,A,B) Y=MATMUL(A,TRANSPOSE(B)) !#define __MATMAT_TT(Y,A,B) Y=TRANSPOSE(MATMUL(B,A))

!#define __PMATMAT_NN(fy,Y,A,B) Y=fyY+MATMUL(A,B) !#define __PMATMAT_TN(fy,Y,A,B) Y=fyY+MATMUL(TRANSPOSE(A),B) !#define __PMATMAT_NT(fy,Y,A,B) Y=fyY+MATMUL(A,TRANSPOSE(B)) !#define __PMATMAT_TT(fy,Y,A,B) Y=fyY+TRANSPOSE(MATMUL(B,A))

!#define __AMATMAT_NN(Y,fa,A,B) Y=faMATMUL(A,B) !#define __AMATMAT_TN(Y,fa,A,B) Y=faMATMUL(TRANSPOSE(A),B) !#define __AMATMAT_NT(Y,fa,A,B) Y=faMATMUL(A,TRANSPOSE(B)) !#define __AMATMAT_TT(Y,fa,A,B) Y=faTRANSPOSE(MATMUL(B,A))

!#define __PAMATMAT_NN(fy,Y,fa,A,B) Y=fyY+faMATMUL(A,B) !#define __PAMATMAT_TN(fy,Y,fa,A,B) Y=fyY+faMATMUL(TRANSPOSE(A),B) !#define __PAMATMAT_NT(fy,Y,fa,A,B) Y=fyY+faMATMUL(A,TRANSPOSE(B)) !#define __PAMATMAT_TT(fy,Y,fa,A,B) Y=fyY+faTRANSPOSE(MATMUL(B,A))

! GEMM does in general Y = fa A^?B^? + fy Y ! with structure: (m x n) = (m x k) (k x n) ! Y=A B : DGEMM('N','N',m,n,k,fa,Amat ,m, Bmat,k, fy,Y,m) ! Y=A^TB : DGEMM('T','N',m,n,k,fa,Amat ,k, Bmat,k, fy,Y,m) ! Y=A B^T : DGEMM('N','T',m,n,k,fa,Amat ,m, Bmat,n, fy,Y,m) ! Y=A^T*B^T : DGEMM('T','T',m,n,k,fa,Amat ,k, Bmat,n, fy,Y,m)

!#define __GENERICMATMAT_NN(fy,Y,fa,A,B) CALL DGEMM('N','N',SIZE(A,1),SIZE(B,2),SIZE(B,1),fa,A,SIZE(A,1),B,SIZE(B,1),fy,Y,SIZE(A,1)) !#define __GENERICMATMAT_TN(fy,Y,fa,A,B) CALL DGEMM('T','N',SIZE(A,2),SIZE(B,2),SIZE(B,1),fa,A,SIZE(B,1),B,SIZE(B,1),fy,Y,SIZE(A,2)) !#define __GENERICMATMAT_NT(fy,Y,fa,A,B) CALL DGEMM('N','T',SIZE(A,1),SIZE(B,1),SIZE(B,2),fa,A,SIZE(A,1),B,SIZE(B,1),fy,Y,SIZE(A,1)) !#define __GENERICMATMAT_TT(fy,Y,fa,A,B) CALL DGEMM('T','T',SIZE(A,2),SIZE(B,1),SIZE(B,2),fa,A,SIZE(B,2),B,SIZE(B,1),fy,Y,SIZE(A,2))

! SIMPLE INTERFACE FOR DGEMM, specifying nrows/ncols of mat A and nrows/ncols of mat B (for any transpose!) ! GEMM does in general Y = fa A^?B^? + fy Y ! with structure: (m x n) = (m x k) (k x n) ! Y=A B : DGEMM('N','N',m,n,k,fa,Amat ,m, Bmat,k, fy,Y,m) ! Y=A^TB : DGEMM('T','N',m,n,k,fa,Amat ,k, Bmat,k, fy,Y,m) ! Y=A B^T : DGEMM('N','T',m,n,k,fa,Amat ,m, Bmat,n, fy,Y,m) ! Y=A^T*B^T : DGEMM('T','T',m,n,k,fa,Amat ,k, Bmat,n, fy,Y,m)

This file depends on

Help

Files dependent on this one

Help

Source Code

!===================================================================================================================================
! Copyright (c) 2025 GVEC Contributors, Max Planck Institute for Plasma Physics
! License: MIT
!===================================================================================================================================
#include "defines.FPP"

!===================================================================================================================================
!>
!!# Module ** rProfile **
!!
!! Defines a 1-D profile in rho^2 via B-Spline knots and coefficients
!===================================================================================================================================
MODULE MODgvec_rProfile_bspl
  ! MODULES
  USE MODgvec_Globals ,ONLY: wp, abort
  USE MODgvec_rProfile_base, ONLY: c_rProfile
  USE sll_m_bsplines  ,ONLY: sll_s_bsplines_new, sll_c_bsplines
  IMPLICIT NONE

  PUBLIC

  TYPE, EXTENDS(c_rProfile) :: t_rProfile_bspl

    INTEGER               :: n_knots !! number of knots, including repeated edge knots
    !INTEGER               :: n_coefs !! number of B-Spline coefficients, part of abstract type
    INTEGER               :: deg = 0
    REAL(wp), ALLOCATABLE :: knots(:)   !! knot values, includinng edge knots
    !REAL(wp), ALLOCATABLE :: coefs(:)   !! B-Spline coefficients, part of abstract type
    CLASS(sll_c_bsplines),ALLOCATABLE :: bspl !! b-spline class

    CONTAINS
    PROCEDURE :: eval_at_rho2        => bsplProfile_eval_at_rho2
    PROCEDURE :: antiderivative      => bsplProfile_antiderivative
    FINAL :: bsplProfile_free

  END TYPE t_rProfile_bspl

  INTERFACE t_rProfile_bspl
      MODULE PROCEDURE bsplProfile_new
  END INTERFACE t_rProfile_bspl


  CONTAINS

  !===================================================================================================================================
  !> initialize the rProfile of type bspline
  !!
  !===================================================================================================================================
  FUNCTION bsplProfile_new(knots,  coefs) RESULT(sf)
    ! INPUT VARIABLES -------------------------!
    REAL(wp),    INTENT(IN) :: knots(:)  !! knots of the B-Spline with repeated start and end points
    REAL(wp),    INTENT(IN) :: coefs(:)  !! B-Spline coefficients
    ! OUTPUT VARIABLES -------------------------!
    TYPE(t_rProfile_bspl) :: sf !! self
    ! LOCAL VARIABLES -------------------------!
    INTEGER :: n_knots, n_coefs
    ! CODE --------------------------------------------------------------------------------------------------------------------------!
    n_knots=SIZE(knots)
    n_coefs=SIZE(coefs)
    sf%deg   = COUNT((ABS(knots-knots(1)).LE.1E-12))-1 ! multiplicity of the first knot determines the degree
    IF(COUNT((ABS(knots-knots(n_knots)).LE.1E-12)).NE.sf%deg+1) THEN
      CALL abort(__STAMP__, &
                 "The Bspline knot sequence need the same multiplicity at the beginning and end (=degree+1).")
    END IF
    IF (n_coefs .NE. n_knots-sf%deg-1) THEN
      CALL abort(__STAMP__, &
                 "Number of Bspline coeffcients must be number of knots - (degree+1)!")
    END IF
    sf%n_knots = n_knots
    sf%n_coefs = n_coefs
    ALLOCATE(sf%knots(1:n_knots), sf%coefs(1:n_coefs))
    sf%knots = knots
    sf%coefs = coefs
    IF (sf%deg>0) THEN
      CALL sll_s_bsplines_new(sf%bspl, sf%deg, .FALSE., &
                              sf%knots(1),sf%knots(n_knots),&
                              size(sf%knots(sf%deg+1:n_knots-sf%deg))-1 , & ! number of knots handed to the library
                              sf%knots(sf%deg+1:n_knots-sf%deg)) ! remove repeated edge knots
    END IF
  END FUNCTION bsplProfile_new

  !===================================================================================================================================
  !> evaluate the n-th derivative of the bsplProfile at position s
  !!
  !===================================================================================================================================
  FUNCTION bsplProfile_eval_at_rho2( sf, rho2, deriv ) RESULT(profile_prime_value)
    ! INPUT VARIABLES -------------------------!
    CLASS(t_rProfile_bspl), INTENT(IN)  :: sf !! self
    REAL(wp)              , INTENT(IN)  :: rho2 !! evaluation point in the toroidal flux coordinate (rho2=phi/phi_edge= rhopos^2)
    INTEGER , OPTIONAL    , INTENT(IN)  :: deriv !! derivative of bspline(rho^2) in rho^2
    ! OUTPUT VARIABLES -------------------------!
    REAL(wp)                         :: profile_prime_value
    ! LOCAL VARIABLES -------------------------!
    REAL(wp)           :: base_eval(0:sf%deg,0:sf%deg) !! value and derivatives of the (deg+1) B-splines that contribute at s_pos
    INTEGER            :: first_non_zero_bspl !! index offset for the coefficients
    INTEGER             :: deriv_case
    ! CODE --------------------------------------------------------------------------------------------------------------------------!
    IF (PRESENT(deriv)) THEN
      deriv_case = deriv
    ELSE
      deriv_case = 0
    END IF
    IF(deriv_case.LE.sf%deg)THEN
      CALL sf%bspl%eval_basis_and_n_derivs(rho2,deriv_case,base_eval(0:deriv_case,:),first_non_zero_bspl)
      profile_prime_value =  SUM(sf%coefs(first_non_zero_bspl:first_non_zero_bspl+sf%deg)*base_eval(deriv_case,:))
    ELSE
      profile_prime_value = 0.0_wp
    END IF
  END FUNCTION bsplProfile_eval_at_rho2

  !===================================================================================================================================
  !> get the exact spline antiderivative, with respect to rho2
  !! the knotspan is increased by an extra multiplicity on both ends, and the new coefficients are computed as
  !! beta(i) = beta(i-1) + alpha(i)*(t(i+degree+1)-t(i))/(degree+1)
  !! From deBoor, "A practical guide to Splines", p.128
  !!
  !===================================================================================================================================
  FUNCTION bsplProfile_antiderivative(sf) RESULT(antideriv)
    ! INPUT VARIABLES -------------------------!
    CLASS(t_rProfile_bspl), INTENT(IN)  :: sf !! self
    ! OUTPUT VARIABLES -------------------------!
    CLASS(c_rProfile),ALLOCATABLE :: antideriv
    ! LOCAL VARIABLES -------------------------!
    REAL(wp) :: coefs(sf%n_coefs+1), knots(sf%n_knots+2), intermid
    INTEGER :: i, n_coefs, n_knots, deg
    ! CODE --------------------------------------------------------------------------------------------------------------------------!
    coefs = 0.0_wp
    knots = -42.0_wp

    n_coefs = sf%n_coefs+1
    n_knots = sf%n_knots+2
    deg = sf%deg+1
    ! increase multiplicity at the edges
    knots(2:n_knots-1) = sf%knots
    knots(1) = sf%knots(1)
    knots(n_knots) = sf%knots(sf%n_knots)

    DO i=1,sf%n_coefs
      intermid = sf%coefs(i)*(sf%knots(i+deg)-sf%knots(i))/deg
      coefs(i+1) = coefs(i) + intermid
    END DO
    antideriv = t_rProfile_bspl(knots,  coefs)
  END FUNCTION bsplProfile_antiderivative

  !===================================================================================================================================
  !> finalize the type rProfile
  !!
  !===================================================================================================================================
  SUBROUTINE bsplProfile_free(sf)
    ! INPUT VARIABLES -------------------------!
    TYPE(t_rProfile_bspl), INTENT(INOUT) :: sf !! self
    ! CODE --------------------------------------------------------------------------------------------------------------------------!
    IF (ALLOCATED(sf%bspl)) THEN
      CALL sf%bspl%free()
      DEALLOCATE(sf%bspl)
    END IF
    SDEALLOCATE(sf%knots)
    SDEALLOCATE(sf%coefs)
  END SUBROUTINE bsplProfile_free

END MODULE MODgvec_rProfile_bspl

rprofile_bspline.F90 Source File

Contents

Modules

Source Code

This file depends on

Files dependent on this one

Source Code