hmap.F90 – GVEC

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

Source Code

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

!===================================================================================================================================
!>
!!# Module ** HMAP new **
!!
!!
!!
!===================================================================================================================================
MODULE MODgvec_hmap
! MODULES
USE MODgvec_c_hmap,     ONLY: c_hmap       ,c_hmap_auxvar
#ifdef PP_WHICH_HMAP
#  if PP_WHICH_HMAP == 1
USE MODgvec_hmap_RZ,    ONLY: t_hmap_RZ    ,t_hmap_RZ_auxvar
#  elif PP_WHICH_HMAP == 3
USE MODgvec_hmap_cyl,   ONLY: t_hmap_cyl   ,t_hmap_cyl_auxvar
#  elif PP_WHICH_HMAP == 10
USE MODgvec_hmap_knot,  ONLY: t_hmap_knot  ,t_hmap_knot_auxvar
#  elif PP_WHICH_HMAP == 20
USE MODgvec_hmap_frenet,ONLY: t_hmap_frenet,t_hmap_frenet_auxvar
#  elif PP_WHICH_HMAP == 21
USE MODgvec_hmap_axisNB,ONLY: t_hmap_axisNB,t_hmap_axisNB_auxvar
#  endif
#else
USE MODgvec_hmap_RZ,    ONLY: t_hmap_RZ    ,t_hmap_RZ_auxvar
USE MODgvec_hmap_cyl,   ONLY: t_hmap_cyl   ,t_hmap_cyl_auxvar
USE MODgvec_hmap_knot,  ONLY: t_hmap_knot  ,t_hmap_knot_auxvar
USE MODgvec_hmap_frenet,ONLY: t_hmap_frenet,t_hmap_frenet_auxvar
USE MODgvec_hmap_axisNB,ONLY: t_hmap_axisNB,t_hmap_axisNB_auxvar
#endif /*defined(PP_WHICH_HMAP)*/

IMPLICIT NONE
PUBLIC


CONTAINS


!===================================================================================================================================
!> initialize the type hmap, also readin parameters here if necessary
!!
!===================================================================================================================================
SUBROUTINE hmap_new( sf, which_hmap,hmap_in)
! MODULES
USE MODgvec_Globals   , ONLY: abort,enter_subregion,exit_subregion
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  INTEGER       , INTENT(IN   ) :: which_hmap         !! input number of field periods
#ifdef PP_WHICH_HMAP
  TYPE(PP_T_HMAP), INTENT(IN),OPTIONAL :: hmap_in       !! if present, copy this hmap
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  TYPE(PP_T_HMAP),ALLOCATABLE,INTENT(INOUT) :: sf !! self
#else
  CLASS(c_hmap), INTENT(IN),OPTIONAL :: hmap_in       !! if present, copy this hmap
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  CLASS(c_hmap),ALLOCATABLE,INTENT(INOUT) :: sf !! self
#endif /*defined(PP_WHICH_HMAP)*/
!===================================================================================================================================
  CALL enter_subregion("hmap")
  IF(.NOT. PRESENT(hmap_in))THEN
    SELECT CASE(which_hmap)
#ifdef PP_WHICH_HMAP
    CASE(PP_WHICH_HMAP)
      sf=PP_T_HMAP()
    CASE DEFAULT
      CALL abort(__STAMP__, &
           "FIXED HMAP TO PP_WHICH_HMAP AT COMPILE TIME,  hmap choice is therefore not compatible  !")
#else
    CASE(1)
      sf=t_hmap_RZ()
    !CASE(2)
    !  ALLOCATE(t_hmap_RphiZ :: sf)
    CASE(3)
      sf=t_hmap_cyl()
    CASE(10)
      sf=t_hmap_knot()
    CASE(20)
      sf=t_hmap_frenet()
    CASE(21)
      sf=t_hmap_axisNB()
    CASE DEFAULT
      CALL abort(__STAMP__, &
           "this hmap choice does not exist  !")
#endif /*defined(PP_WHICH_HAP)*/
    END SELECT
    sf%which_hmap=which_hmap
  ELSE
    IF(which_hmap.NE.hmap_in%which_hmap) CALL abort(__STAMP__, &
       "hmap_in does not coincide with requested hmap in hmap_new")
    ALLOCATE(sf,source=hmap_in)
  END IF
  CALL exit_subregion("hmap")
END SUBROUTINE hmap_new

!===================================================================================================================================
!> initialize the  hmap auxiliary variables, depends on hmap type
!!
!===================================================================================================================================
SUBROUTINE hmap_new_auxvar(hmap,zeta,xv,do_2nd_der)
! MODULES
USE MODgvec_Globals   , ONLY: abort,wp
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
#ifdef PP_WHICH_HMAP
  TYPE(PP_T_HMAP),INTENT(IN) :: hmap
#else
  CLASS(c_hmap),  INTENT(IN) :: hmap
#endif
  REAL(wp)     ,  INTENT(IN) :: zeta(:)
  LOGICAL      ,  INTENT(IN) :: do_2nd_der
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
#ifdef PP_WHICH_HMAP
  TYPE(PP_T_HMAP_AUXVAR),ALLOCATABLE,INTENT(INOUT) :: xv(:) !! self
#else
  CLASS(c_hmap_auxvar),ALLOCATABLE,INTENT(INOUT) :: xv(:) !! self
#endif
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER :: i,nzeta
!===================================================================================================================================
  nzeta=SIZE(zeta)
#ifdef PP_WHICH_HMAP
  IF(hmap%which_hmap .NE. PP_WHICH_HMAP) CALL abort(__STAMP__, &
           "FIXED HMAP TO PP_WHICH_HMAP AT COMPILE TIME,  which_hmap choice is therefore not compatible  !")
  ALLOCATE(PP_T_HMAP_AUXVAR :: xv(nzeta))
  !$OMP PARALLEL DO &
  !$OMP   SCHEDULE(STATIC) DEFAULT(SHARED) PRIVATE(i)
  DO i=1,nzeta
    xv(i)= PP_T_HMAP_AUXVAR(hmap,zeta(i),do_2nd_der)
  END DO !i
  !$OMP END PARALLEL DO

#else

  SELECT TYPE(hmap)
  CLASS IS(t_hmap_RZ)
    ALLOCATE(t_hmap_RZ_auxvar :: xv(nzeta))
    SELECT TYPE(xv)
    TYPE IS(t_hmap_RZ_auxvar)
      !$OMP PARALLEL DO &
      !$OMP   SCHEDULE(STATIC) DEFAULT(SHARED) PRIVATE(i)
      DO i=1,nzeta
        xv(i)= t_hmap_RZ_auxvar(hmap,zeta(i),do_2nd_der)
      END DO !i
      !$OMP END PARALLEL DO
    END SELECT !TYPE(xv)
  CLASS IS(t_hmap_cyl)
    ALLOCATE(t_hmap_cyl_auxvar :: xv(nzeta))
    SELECT TYPE(xv)
    TYPE IS(t_hmap_cyl_auxvar)
      !$OMP PARALLEL DO &
      !$OMP   SCHEDULE(STATIC) DEFAULT(SHARED) PRIVATE(i)
      DO i=1,nzeta
        xv(i)= t_hmap_cyl_auxvar(hmap,zeta(i),do_2nd_der)
      END DO !i
      !$OMP END PARALLEL DO
    END SELECT !TYPE(xv)
  CLASS IS(t_hmap_knot)
    ALLOCATE(t_hmap_knot_auxvar :: xv(nzeta))
    SELECT TYPE(xv)
    TYPE IS(t_hmap_knot_auxvar)
      !$OMP PARALLEL DO &
      !$OMP   SCHEDULE(STATIC) DEFAULT(SHARED) PRIVATE(i)
      DO i=1,nzeta
        xv(i)= t_hmap_knot_auxvar(hmap,zeta(i),do_2nd_der)
      END DO !i
      !$OMP END PARALLEL DO
    END SELECT !TYPE(xv)
  CLASS IS(t_hmap_frenet)
    ALLOCATE(t_hmap_frenet_auxvar :: xv(nzeta))
    SELECT TYPE(xv)
    TYPE IS(t_hmap_frenet_auxvar)
      !$OMP PARALLEL DO &
      !$OMP   SCHEDULE(STATIC) DEFAULT(SHARED) PRIVATE(i)
      DO i=1,nzeta
        xv(i)= t_hmap_frenet_auxvar(hmap,zeta(i),do_2nd_der)
      END DO !i
      !$OMP END PARALLEL DO
    END SELECT !TYPE(xv)
  CLASS IS(t_hmap_axisNB)
    ALLOCATE(t_hmap_axisNB_auxvar :: xv(nzeta))
    SELECT TYPE(xv)
    TYPE IS(t_hmap_axisNB_auxvar)
      !$OMP PARALLEL DO &
      !$OMP   SCHEDULE(STATIC) DEFAULT(SHARED) PRIVATE(i)
      DO i=1,nzeta
        xv(i)= t_hmap_axisNB_auxvar(hmap,zeta(i),do_2nd_der)
      END DO !i
      !$OMP END PARALLEL DO
    END SELECT !TYPE(xv)
  CLASS DEFAULT
    CALL abort(__STAMP__, &
          "hmap_new_auxvar: this hmap class is not implemented  !")
  END SELECT
#endif /*defined(PP_WHICH_HAP)*/

END SUBROUTINE hmap_new_auxvar

END MODULE MODgvec_hmap

hmap.F90 Source File

Contents

Modules

Source Code

This file depends on

Files dependent on this one

Source Code