rprofile_polynomial.F90 Source File

!!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

sourcefile~~rprofile_polynomial.f90~~EfferentGraph sourcefile~rprofile_polynomial.f90 rprofile_polynomial.F90 sourcefile~c_rprofile.f90 c_rprofile.F90 sourcefile~rprofile_polynomial.f90->sourcefile~c_rprofile.f90 sourcefile~globals.f90 globals.F90 sourcefile~rprofile_polynomial.f90->sourcefile~globals.f90 sourcefile~c_rprofile.f90->sourcefile~globals.f90

Files dependent on this one

sourcefile~~rprofile_polynomial.f90~~AfferentGraph sourcefile~rprofile_polynomial.f90 rprofile_polynomial.F90 sourcefile~mhd3d.f90 mhd3d.F90 sourcefile~mhd3d.f90->sourcefile~rprofile_polynomial.f90 sourcefile~gvec_post.f90 gvec_post.F90 sourcefile~gvec_post.f90->sourcefile~mhd3d.f90 sourcefile~rungvec.f90 rungvec.F90 sourcefile~rungvec.f90->sourcefile~mhd3d.f90 sourcefile~state.f90 state.F90 sourcefile~state.f90->sourcefile~mhd3d.f90 sourcefile~gvec.f90 gvec.F90 sourcefile~gvec.f90->sourcefile~rungvec.f90 sourcefile~run.f90 run.F90 sourcefile~run.f90->sourcefile~rungvec.f90

Source Code

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

!===================================================================================================================================
!>
!!# Module ** polyProfile **
!!
!! Defines a 1-D profile in rho^2 via a power polynomial.
!===================================================================================================================================
MODULE MODgvec_rProfile_poly
  ! MODULES
  USE MODgvec_Globals ,ONLY: wp
  USE MODgvec_rProfile_base, ONLY: c_rProfile, poly_derivative_prefactor
  IMPLICIT NONE

  PUBLIC

  TYPE, EXTENDS(c_rProfile) :: t_rProfile_poly
    !INTEGER               :: n_coefs !! number of polynomial coefficients, part of abstract type
    INTEGER               :: deg = 0
    !REAL(wp), ALLOCATABLE :: coefs(:)   !! polynomial coefficients, part of abstract type
    CONTAINS
    PROCEDURE :: eval_at_rho2        => polyProfile_eval_at_rho2
    PROCEDURE :: antiderivative      => polyProfile_antiderivative
    FINAL     :: polyProfile_free

  END TYPE t_rProfile_poly

  INTERFACE t_rProfile_poly
      MODULE PROCEDURE polyProfile_new
  END INTERFACE t_rProfile_poly

  CONTAINS

  FUNCTION polyProfile_new(coefs) RESULT(sf)
    ! INPUT VARIABLES -------------------------!
    REAL,    INTENT(IN) :: coefs(:)  !! B-Spline coefficients
    ! OUTPUT VARIABLES -------------------------!
    TYPE(t_rProfile_poly) :: sf
    ! LOCAL VARIABLES -------------------------!
    INTEGER :: n_coefs
    ! CODE --------------------------------------------------------------------------------------------------------------------------!
    n_coefs=SIZE(coefs)
    sf%deg   = n_coefs-1
    sf%n_coefs = n_coefs
    ALLOCATE(sf%coefs(1:n_coefs))
    sf%coefs = coefs
  END FUNCTION polyProfile_new

  !===================================================================================================================================
  !> evaluate the n-th derivative of a power polynomial
  !!
  !===================================================================================================================================
  FUNCTION polyProfile_eval_at_rho2(sf, rho2, deriv) RESULT(profile_prime_value)
    ! MODULES
    USE MODgvec_Globals, ONLY: Eval1DPoly,Eval1DPoly_deriv
    ! INPUT VARIABLES -------------------------!
    CLASS(t_rProfile_poly), INTENT(IN)  :: sf !! self
    REAL(wp)              , INTENT(IN)  :: rho2 !! evaluation point in the toroidal flux coordinate (rho2=phi/phi_edge= spos^2)
    INTEGER , OPTIONAL    , INTENT(IN)  :: deriv
    ! OUTPUT VARIABLES -------------------------!
    REAL(wp)                         :: profile_prime_value
    ! LOCAL VARIABLES -------------------------!
    REAL(wp)                         :: prefactor
    INTEGER                          :: d
    INTEGER                          :: deriv_case
    ! CODE --------------------------------------------------------------------------------------------------------------------------!
    IF (PRESENT(deriv)) THEN
        deriv_case = deriv
    ELSE
        deriv_case = 0
    END IF

    IF (deriv_case>sf%deg) THEN
        profile_prime_value = 0.0_wp
    ELSE IF (deriv_case==0) THEN
        profile_prime_value = EVAL1DPOLY(sf%n_coefs, sf%coefs, rho2)
    ELSE IF (deriv_case==1) THEN
        profile_prime_value = EVAL1DPOLY_deriv(sf%n_coefs, sf%coefs, rho2)
    ELSE
        profile_prime_value = 0.0_wp
        DO d=sf%deg+1, deriv_case+1,-1
            prefactor=poly_derivative_prefactor(d-1,deriv_case)
            profile_prime_value = profile_prime_value*rho2+prefactor*sf%coefs(d)
        END DO
    END IF
  END FUNCTION polyProfile_eval_at_rho2

  !===================================================================================================================================
  !> get the exact polynomial antiderivative, with respect to rho2
  !!
  !===================================================================================================================================
  FUNCTION polyProfile_antiderivative(sf) RESULT(antideriv)
    ! INPUT VARIABLES -------------------------!
    CLASS(t_rProfile_poly), INTENT(IN)  :: sf !! self
    ! OUTPUT VARIABLES -------------------------!
    CLASS(c_rProfile),ALLOCATABLE :: antideriv
    ! LOCAL VARIABLES -------------------------!
    REAL(wp) :: coefs(sf%n_coefs+1)
    INTEGER :: i
    ! CODE --------------------------------------------------------------------------------------------------------------------------!
    coefs = 0.0_wp
    DO i=1,sf%n_coefs
      coefs(i+1) = sf%coefs(i)/REAL(i,wp)
    END DO
    antideriv = t_rProfile_poly(coefs)
  END FUNCTION polyProfile_antiderivative

  !===================================================================================================================================
  !> finalize the type rProfile
  !!
  !===================================================================================================================================
  SUBROUTINE polyProfile_free(sf)
    ! INPUT VARIABLES -------------------------!
    TYPE(t_rProfile_poly), INTENT(INOUT)  :: sf !! self
    ! CODE --------------------------------------------------------------------------------------------------------------------------!
    SDEALLOCATE(sf%coefs)
  END SUBROUTINE polyProfile_free

END MODULE MODgvec_rProfile_poly