c_rprofile.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 **
!!
!! Abstract class for radial profiles
!===================================================================================================================================
MODULE MODgvec_rProfile_base
  ! MODULES
  USE MODgvec_Globals ,ONLY: wp, abort
  IMPLICIT NONE

  PUBLIC

  TYPE, ABSTRACT :: c_rProfile
      INTEGER               :: n_coefs
      REAL(wp), ALLOCATABLE :: coefs(:)

      contains

      PROCEDURE(i_fun_eval_at_rho2), DEFERRED :: eval_at_rho2
      PROCEDURE(i_fun_antiderivative), DEFERRED :: antiderivative

      PROCEDURE :: eval_at_rho => rProfile_eval_at_rho
      ! hard coded derivatives with respect to rho=sqrt(phi/phi_edge)
      PROCEDURE :: rProfile_drho2
      PROCEDURE :: rProfile_drho3
      PROCEDURE :: rProfile_drho4

  end type c_rProfile

  ABSTRACT INTERFACE

      FUNCTION i_fun_eval_at_rho2( sf, rho2, deriv ) RESULT(profile_value)
          IMPORT c_rProfile
          IMPORT wp
          CLASS(c_rProfile), INTENT(IN)  :: sf
          REAL(wp)         , INTENT(IN)  :: rho2
          INTEGER, OPTIONAL, INTENT(IN)  :: deriv
          REAL(wp)                       :: profile_value
      END FUNCTION i_fun_eval_at_rho2

      FUNCTION i_fun_antiderivative(sf) RESULT(antideriv)
          IMPORT c_rProfile
          IMPORT wp
          CLASS(c_rProfile), INTENT(IN)   :: sf
          CLASS(c_rProfile), ALLOCATABLE  :: antideriv
      END FUNCTION i_fun_antiderivative
  END INTERFACE

  CONTAINS
  !===================================================================================================================================
  !> calculate the prefactor for the d-th coefficient of the n-th derivative of a polynomial
  !!
  !===================================================================================================================================
  PURE FUNCTION poly_derivative_prefactor(D,deriv) RESULT(prefactor)
    INTEGER, INTENT(IN) :: D,deriv
    INTEGER :: i
    REAL(wp) :: prefactor
    ! CODE --------------------------------------------------------------------------------------------------------------------------!
    prefactor = 1.0_wp
    DO i=D-deriv+1,D
        prefactor = prefactor*i
    END DO
  END FUNCTION poly_derivative_prefactor

  !===================================================================================================================================
  !> evaluate the n-th derivative of (rho^2) with respect to rho ~sqrt(magnetic flux).
  !!
  !===================================================================================================================================
  PURE FUNCTION rho2_derivative(rho,deriv) RESULT(rho2_prime)
    ! INPUT VARIABLES -------------------------!
    REAL(wp), INTENT(IN) :: rho                !! rho position rho ~sqrt(magnetic flux)
    INTEGER, INTENT(IN)  :: deriv  !! derivative in rho
    ! OUTPUT VARIABLES -------------------------!
    REAL(wp) :: rho2_prime     !!n-th derivative of rho2 with respect to rho ~sqrt(magnetic flux).
    ! CODE --------------------------------------------------------------------------------------------------------------------------!
    IF (deriv>2) THEN
      rho2_prime = 0.0_wp
  ELSE
      rho2_prime = poly_derivative_prefactor(2,deriv)*rho**(2-deriv)
  END IF
  END FUNCTION rho2_derivative

  !===================================================================================================================================
  !> evaluate the 2nd derivative of a radial profile with respect to rho ~sqrt(magnetic flux).
  !!
  !===================================================================================================================================
  FUNCTION rProfile_drho2(sf, rho) RESULT(derivative)
    ! INPUT VARIABLES -------------------------!
    CLASS(c_rProfile)  :: sf !! self
    REAL(wp), INTENT(IN) :: rho                !! rho position rho ~sqrt(magnetic flux)
    ! OUTPUT VARIABLES -------------------------!
    REAL(wp) :: derivative   !! 2nd derivative of a radial profile with respect to rho ~sqrt(magnetic flux).
    ! LOCAL VARIABLES -------------------------!
    REAL(wp) :: rho2
    ! CODE --------------------------------------------------------------------------------------------------------------------------!
    rho2 = rho2_derivative(rho,deriv=0)
    ! d^2/dx^2 f(g(x)) = [g'(x)]^2 * f''(g(x))+g''(x) * f'(g(x))
    derivative =  rho2_derivative(rho,deriv=1)**2*sf%eval_at_rho2(rho2, deriv=2) &
                + rho2_derivative(rho,deriv=2)*   sf%eval_at_rho2(rho2, deriv=1)
  END FUNCTION rProfile_drho2

  !===================================================================================================================================
  !> evaluate the 3rd derivative of a radial profile with respect to rho ~sqrt(magnetic flux).
  !!
  !===================================================================================================================================
  FUNCTION rProfile_drho3(sf, rho) RESULT(derivative)
    ! INPUT VARIABLES -------------------------!
    CLASS(c_rProfile)  :: sf !! self
    REAL(wp), INTENT(IN) :: rho                !! rho position rho ~sqrt(magnetic flux)
    ! OUTPUT VARIABLES -------------------------!
    REAL(wp) :: derivative          !! 3rd derivative of a radial profile with respect to rho ~sqrt(magnetic flux)
    ! LOCAL VARIABLES -------------------------!
    REAL(wp) :: rho2
    ! CODE --------------------------------------------------------------------------------------------------------------------------!
    rho2 = rho2_derivative(rho,deriv=0)
    ! d^3/dx^3 f(g(x)) = 3g'(x) * f''(g(x)) * g''(x) + [g'(x)]^3*f'''(g(x)) + f'(x) * g'''(x)
    derivative = 3*rho2_derivative(rho,deriv=1)*   sf%eval_at_rho2(rho2, deriv=2)*rho2_derivative(rho,deriv=2) &
                +  rho2_derivative(rho,deriv=1)**3*sf%eval_at_rho2(rho2, deriv=3) &
                +  rho2_derivative(rho,deriv=3)*   sf%eval_at_rho2(rho2, deriv=1)
  END FUNCTION rProfile_drho3

  !===================================================================================================================================
  !> evaluate the 4th derivative of a radial profile with respect to rho ~sqrt(magnetic flux)
  !!
  !===================================================================================================================================
  FUNCTION rProfile_drho4(sf, rho) RESULT(derivative)
    ! INPUT VARIABLES -------------------------!
    CLASS(c_rProfile) :: sf !! self
    REAL(wp), INTENT(IN) :: rho                !! rho position rho ~sqrt(magnetic flux)
    ! OUTPUT VARIABLES -------------------------!
    REAL(wp) :: derivative   !! 4th derivative of a radial profile with respect to rho ~sqrt(magnetic flux).
    ! LOCAL VARIABLES -------------------------!
    REAL(wp) :: rho2
    ! CODE --------------------------------------------------------------------------------------------------------------------------!
    rho2 = rho2_derivative(rho,deriv=0)
    ! d^4/dx^4 f(g(x)) = f''''(g(x))g'(x)**4
    !                   + 6f'''(g(x))g''(x)g'(x)^2
    !                   + 3f''(g(x))g''(x)^2
    !                   + 4f''(g(x))g'''(x)g'(x)
    !                   + f'(g(x))g''''(x)
    derivative =    sf%eval_at_rho2(rho2, deriv=4)*rho2_derivative(rho,deriv=1)**4 &
                + 6*sf%eval_at_rho2(rho2, deriv=3)*rho2_derivative(rho,deriv=2)*rho2_derivative(rho,deriv=1)**2 &
                + 3*sf%eval_at_rho2(rho2, deriv=2)*rho2_derivative(rho,deriv=2)**2 &
                + 4*sf%eval_at_rho2(rho2, deriv=2)*rho2_derivative(rho,deriv=3)*rho2_derivative(rho,deriv=1) &
                +   sf%eval_at_rho2(rho2, deriv=1)*rho2_derivative(rho,deriv=4)
  END FUNCTION rProfile_drho4

  !===================================================================================================================================
  !> evaluate the n-th derivative of a radial profile with respect to rho ~sqrt(magnetic flux).
  !! NOTE: n has to be in [0,4] due to an explicit implementation of the product rule.
  !===================================================================================================================================
  FUNCTION rProfile_eval_at_rho(sf, rho, deriv) RESULT(derivative)
    ! INPUT VARIABLES -------------------------!
    CLASS(c_rProfile)  :: sf !! self
    REAL(wp), INTENT(IN) :: rho                !! rho position rho ~sqrt(magnetic flux)
    INTEGER,  OPTIONAL, INTENT(IN)   :: deriv  !! derivative in rho
    ! OUTPUT VARIABLES -------------------------!
    REAL(wp) :: derivative  !! n-th derivative of a radial profile with respect to rho ~sqrt(magnetic flux).
    ! LOCAL VARIABLES -------------------------!
    INTEGER :: deriv_case
    REAL(wp) :: rho2
    ! CODE --------------------------------------------------------------------------------------------------------------------------!
    IF (PRESENT(deriv)) THEN
        deriv_case = deriv
    ELSE
        deriv_case = 0
    END IF

    rho2 = rho2_derivative(rho,deriv=0)
    SELECT CASE(deriv_case)
    CASE(0)
        derivative = sf%eval_at_rho2(rho2, deriv=0)
    CASE(1)
        derivative = sf%eval_at_rho2(rho2,deriv=1)*rho2_derivative(rho,deriv=1)
    CASE(2)
        derivative = sf%rProfile_drho2(rho)
    CASE(3)
        derivative = sf%rProfile_drho3(rho)
    CASE(4)
        derivative = sf%rProfile_drho4(rho)
    CASE DEFAULT
        CALL abort(__STAMP__,&
            "error in rprofile: derivatives higher than 4 with respect to rho=sqrt(phi/phi_edge) are not implemented!")
    END SELECT
  END FUNCTION rProfile_eval_at_rho

END MODULE MODgvec_rProfile_base

c_rprofile.F90 Source File

Contents

Modules

Source Code

This file depends on

Files dependent on this one

Source Code