!!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)
!=================================================================================================================================== ! Copyright (c) 2025 GVEC Contributors, Max Planck Institute for Plasma Physics ! License: MIT !=================================================================================================================================== #include "defines.FPP" !=================================================================================================================================== !> !!# Module **Read in the boundary from a specifically defined NETCDF file ** !! !! !=================================================================================================================================== MODULE MODgvec_boundaryFromFile ! MODULES USE MODgvec_Globals, ONLY:wp,abort,UNIT_stdOut USE MODgvec_io_netcdf, ONLY:t_ncfile IMPLICIT NONE PRIVATE TYPE :: t_boundaryFromFile !--------------------------------------------------------------------------------------------------------------------------------- LOGICAL :: initialized=.FALSE. !! set to true in init, set to false in free !--------------------------------------------------------------------------------------------------------------------------------- INTEGER :: nfp !! number of field periods INTEGER :: m_max !! maximum number of fourier modes of the boundary in theta direction INTEGER :: n_max !! maximum number of fourier modes of the boundary in zeta direction (one field period only!) INTEGER :: ntheta,nzeta !! number of interpolation points in theta and zeta (one field period, half grid!!) INTEGER :: lasym !! =0: symmetric, =1: asymmetric fourier series REAL(wp),ALLOCATABLE :: theta(:) !! theta positions [0,2pi), should be half grid! REAL(wp),ALLOCATABLE :: zeta(:) !! zeta positions [0,2pi/nfp) should be on half grid! REAL(wp),ALLOCATABLE :: X(:,:) ,Y(:,:)!! boundary data X/Y positions X[i, j]=X(theta[i],zeta[j]), !! Y[i, j]=Y(theta[i],zeta[j]), i=0...ntheta-1,j=0...nzeta-1 CLASS(t_ncfile),ALLOCATABLE :: nc !! container for netcdf-file CHARACTER(LEN=1024) :: ncfile=" " !! name of netcdf file with axis information CONTAINS PROCEDURE :: init => bff_init PROCEDURE :: convert_to_modes => bff_convert_to_modes PROCEDURE :: free => bff_free END TYPE t_boundaryFromFile INTERFACE boundaryFromFile_new MODULE PROCEDURE boundaryFromFile_new END INTERFACE PUBLIC::t_boundaryFromFile,boundaryFromFile_new !=================================================================================================================================== CONTAINS !=================================================================================================================================== !> Allocate class and call init !! !=================================================================================================================================== SUBROUTINE boundaryFromFile_new( sf,fileString) ! MODULES IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES CHARACTER(LEN=*) , INTENT(IN ) :: fileString !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES CLASS(t_boundaryFromFile), ALLOCATABLE,INTENT(INOUT) :: sf !! self !=================================================================================================================================== ALLOCATE(t_boundaryFromFile :: sf) CALL sf%init(fileString) END SUBROUTINE boundaryFromFile_new !=================================================================================================================================== !> initialize class: read file and save data to class structure !! !=================================================================================================================================== SUBROUTINE bff_init(sf,fileString) ! MODULES USE MODgvec_io_netcdf, ONLY:ncfile_init IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES CHARACTER(LEN=*) , INTENT(IN ) :: fileString !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES CLASS(t_boundaryFromFile), INTENT(INOUT) :: sf !! self !=================================================================================================================================== WRITE(UNIT_stdOut,'(A)')' READ BOUNDARY FROM NETCDF FILE "'//TRIM(FileString)//'" ...' sf%ncfile=TRIM(FileString) CALL ncfile_init(sf%nc,sf%ncfile,"r") CALL READNETCDF(sf) sf%initialized=.TRUE. END SUBROUTINE bff_init !=================================================================================================================================== !> READ FROM SPECIFIC NETCDF FILE: general data and "boundary" group !! ======= HEADER OF THE NETCDF FILE VERSION 3.1 =================================================================================== !! === FILE DESCRIPTION: !! * axis, normal and binormal of the frame are given in cartesian coordinates along the curve parameter zeta [0,2pi]. !! * The curve is allowed to have a field periodicity NFP, but the curve must be provided on a full turn. !! * The data is given in REAL SPACE, sampled along equidistant zeta point positions: !! zeta(i)=(i+0.5)/nzeta * (2pi/NFP), i=0,...,nzeta-1 !! always shifted by (2pi/NFP) for the next field period. !! Thus the number of points along the axis for a full turn is NFP*nzeta !! * definition of the axis-following frame in cartesian coordinates ( boundary surface at rho=1): !! !! {x,y,z}(rho,theta,zeta)={x,y,z}(zeta) + X(rho,theta,zeta)*N_{x,y,z}(zeta)+Y(rho,theta,zeta)*B_{x,y,z}(zeta) !! !! === DATA DESCRIPTION !! - general data !! * NFP: number of field periods !! * VERSION: version number as integer: V3.0 => 300 !! - axis/ data group: !! * 'axis/n_max' : maximum mode number in zeta (in one field period) !! * 'axis/nzeta' : number of points along the axis, in one field period (>=2*n_max+1) !! * 'axis/zeta(:)' : zeta positions, 1D array of size 'axis/nzeta', for one field period. zeta[i]=zeta[1] + (i-1)/nzeta*(2pi/nfp), i=1,..nzeta, zeta[1] is arbitrary !! * 'axis/xyz(::)' : cartesian positions along the axis for ONE FULL TURN, 2D array of size (3,NFP* nzeta ), sampled at zeta positions, must exclude the endpoint !! xyz[:,j+fp*nzeta]=axis(zeta[j]+fp*2pi/NFP), for j=0,..nzeta-1 and fp=0,...,NFP-1 !! * 'axis/Nxyz(::)': cartesian components of the normal vector of the axis frame, 2D array of size (3, NFP* nzeta), evaluated analogously to the axis !! * 'axis/Bxyz(::)': cartesian components of the bi-normal vector of the axis frame, 2D array of size (3, NFP*nzeta), evaluated analogously to the axis !! - boundary data group: !! * 'boundary/m_max' : maximum mode number in theta !! * 'boundary/n_max' : maximum mode number in zeta (in one field period) !! * 'boundary/lasym' : asymmetry, logical. !! if lasym=0, boundary surface position X,Y in the N-B plane of the axis frame can be represented only with !! X(theta,zeta)=sum X_mn*cos(m*theta-n*NFP*zeta), with {m=0,n=0...n_max},{m=1...m_max,n=-n_max...n_max} !! Y(theta,zeta)=sum Y_mn*sin(m*theta-n*NFP*zeta), with {m=0,n=1...n_max},{m=1...m_max,n=-n_max...n_max} !! if lasym=1, full fourier series is taken for X,Y !! * 'boundary/ntheta' : number of points in theta (>=2*m_max+1) !! * 'boundary/nzeta' : number of points in zeta (>=2*n_max+1), can be different to 'axis/nzeta' !! * 'boundary/theta(:)' : theta positions, 1D array of size 'boundary_ntheta', on half grid! theta(i)=(i+0.5)/ntheta*(2pi), i=0,...ntheta-1 !! * 'boundary/zeta(:)' : zeta positions, 1D array of size 'boundary/nzeta', for one field period! zeta[i]=zeta[1] + (i-1)/nzeta*(2pi/nfp), i=1,..nzeta, zeta[1] is arbitrary !! * 'boundary/X(::)', !! 'boundary/Y(::)' : boundary position X,Y in the N-B plane of the axis frame, in one field period, 2D array of size(ntheta, nzeta), with !! X[i, j]=X(theta[i],zeta[j]) !! Y[i, j]=Y(theta[i],zeta[j]), i=0...ntheta-1,j=0...nzeta-1 !=================================================================================================================================== SUBROUTINE ReadNETCDF(sf) USE MODgvec_io_netcdf IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT/OUTPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES CLASS(t_boundaryFromFile), INTENT(INOUT) :: sf !! self !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES !=================================================================================================================================== CALL sf%nc%get_scalar("NFP",intout=sf%nfp) CALL sf%nc%get_scalar("boundary/lasym",intout=sf%lasym) CALL sf%nc%get_scalar("boundary/ntheta",intout=sf%ntheta) CALL sf%nc%get_scalar("boundary/nzeta",intout=sf%nzeta) IF(sf%nc%var_exists("boundary/m_max"))THEN CALL sf%nc%get_scalar("boundary/m_max",intout=sf%m_max) sf%m_max=MIN(sf%m_max,(sf%ntheta-1)/2) !maximum mode number based on number of interpolation points ntheta>=2*m_max+1 ELSE sf%m_max=(sf%ntheta-1)/2 !maximum mode number based on number of interpolation points ntheta>=2*m_max+1 WRITE(UNIT_stdOut,'(6X,A,I8)')'"boundary/m_max" not found, set to: ',sf%m_max END IF IF(sf%nc%var_exists("boundary/n_max"))THEN CALL sf%nc%get_scalar("boundary/n_max",intout=sf%n_max) sf%n_max = MIN(sf%n_max,(sf%nzeta-1)/2) !maximum mode number based on number of interpolation points nzeta>=2*n_max+1 ELSE sf%n_max=(sf%nzeta-1)/2 !maximum mode number based on number of interpolation points nzeta>=2*n_max+1 END IF WRITE(UNIT_stdOut,'(6X,A,2I4)')'" boundary (m_max,n_max)" is set to: ',sf%m_max,sf%n_max ALLOCATE(sf%theta(sf%ntheta)) CALL sf%nc%get_array("boundary/theta(:)",realout_1d=sf%theta) ALLOCATE(sf%zeta(sf%nzeta)) CALL sf%nc%get_array("boundary/zeta(:)",realout_1d=sf%zeta) ALLOCATE(sf%X(sf%ntheta,sf%nzeta)) CALL sf%nc%get_array("boundary/X(::)",realout_2d=sf%X) ALLOCATE(sf%Y(sf%ntheta,sf%nzeta)) CALL sf%nc%get_array("boundary/Y(::)",realout_2d=sf%Y) CALL sf%nc%closefile() END SUBROUTINE READNETCDF !=================================================================================================================================== !> convert from interpolation points X=> X1_b, Y=> X2_b to fourier modes, given from the input fbase !! convert to maximum allowable number of modes (ntheta>=2*m_max+1, nzeta>=2*n_max+1) !! the final m_max/n_max can be smaller or larger. If larger, a change of base is necessary !! !=================================================================================================================================== SUBROUTINE bff_convert_to_modes(sf,x1_fbase_in,x2_fbase_in,X1_b,X2_b,scale_minor_radius) ! MODULES USE MODgvec_fbase ,ONLY: t_fbase, sin_cos_map IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES TYPE(t_fbase), INTENT(IN):: x1_fbase_in,x2_fbase_in REAL(wp), INTENT(IN) :: scale_minor_radius !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES REAL(wp), INTENT(INOUT) :: X1_b(x1_fbase_in%modes),X2_b(x2_fbase_in%modes) CLASS(t_boundaryFromFile), INTENT(INOUT) :: sf !! self !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES TYPE(t_fBase) :: X_fbase,Y_fbase INTEGER :: i,nIP,mIP,mn_max_pts(2) REAL(wp) :: xn(2,sf%ntheta*sf%nzeta) REAL(wp),ALLOCATABLE :: xydofs(:,:),X12dofs(:,:) !=================================================================================================================================== WRITE(UNIT_stdOut,'(A)')' CONVERT BOUNDARY FROM POINTS TO MODES:' IF(sf%nfp.NE.x1_fbase_in%nfp)CALL abort(__STAMP__,& " error in convert boundary to modes, nfp from GVEC parameterfile does not match to nfp of boundary file",& TypeInfo="InvalidParameterError") mn_max_pts(1:2)=(/(sf%ntheta-1)/2,(sf%nzeta-1)/2/) WRITE(UNIT_stdOut,'(6X,2(A,I4),2A)')' from X,Y(ntheta= ',sf%ntheta,', nzeta= ', & sf%nzeta, '), lasym=',MERGE("symmetric","asymetric",(sf%lasym.EQ.0)) WRITE(UNIT_stdOut,'(6x,2(3A,2(2I4,A)))') ' => X1 ',sin_cos_map(x1_fbase_in%sin_cos), & ', (m_max,n_max)= (',mn_max_pts,')=>(',x1_fbase_in%mn_max,')', & ' , X2 ',sin_cos_map(x2_fbase_in%sin_cos), & ', (m_max,n_max)= (',mn_max_pts,')=>(',x2_fbase_in%mn_max,')' IF(ALL(x1_fbase_in%mn_max.LE.mn_max_pts))THEN !X1_base is smaller/equal X_fbase = t_fBase(x1_fbase_in%mn_max, (/sf%ntheta,sf%nzeta/), & sf%nfp, sin_cos_map(x1_fbase_in%sin_cos), x1_fbase_in%exclude_mn_zero) X1_b = X_fbase%initDOF(RESHAPE(sf%X*scale_minor_radius,(/sf%ntheta*sf%nzeta/)) ,thet_zeta_start=(/sf%theta(1),sf%zeta(1)/)) ELSE X_fbase = t_fBase(mn_max_pts, (/sf%ntheta,sf%nzeta/), & sf%nfp, sin_cos_map(x1_fbase_in%sin_cos), x1_fbase_in%exclude_mn_zero) ALLOCATE(xydofs(1,1:X_fbase%modes),X12dofs(1,1:x1_fbase_in%modes)) xydofs(1,:) = X_fbase%initDOF(RESHAPE(sf%X*scale_minor_radius,(/sf%ntheta*sf%nzeta/)),thet_zeta_start=(/sf%theta(1),sf%zeta(1)/)) CALL x1_fbase_in%change_base(X_fbase,1,xydofs,X12dofs) X1_b=X12dofs(1,:) DEALLOCATE(xydofs,X12dofs) END IF IF(ALL(x2_fbase_in%mn_max.LE.mn_max_pts))THEN !X2_base is smaller/equal Y_fbase = t_fBase(x2_fbase_in%mn_max, (/sf%ntheta,sf%nzeta/), & sf%nfp, sin_cos_map(x2_fbase_in%sin_cos), x2_fbase_in%exclude_mn_zero) X2_b = Y_fbase%initDOF(RESHAPE(sf%Y*scale_minor_radius,(/sf%ntheta*sf%nzeta/)) ,thet_zeta_start=(/sf%theta(1),sf%zeta(1)/)) ELSE Y_fbase = t_fBase(mn_max_pts, (/sf%ntheta,sf%nzeta/), & sf%nfp, sin_cos_map(x2_fbase_in%sin_cos), x2_fbase_in%exclude_mn_zero) ALLOCATE(xydofs(1,1:Y_fbase%modes),X12dofs(1,1:x2_fbase_in%modes)) xydofs(1,:) = Y_fbase%initDOF(RESHAPE(sf%Y*scale_minor_radius,(/sf%ntheta*sf%nzeta/)),thet_zeta_start=(/sf%theta(1),sf%zeta(1)/)) CALL x2_fbase_in%change_base(Y_fbase,1,xydofs,X12dofs) X2_b=X12dofs(1,:) END IF !evaluate at interpolation points and check the error i=0 DO nIP=1,sf%nzeta DO mIP=1,sf%ntheta i=i+1 xn(1,i)=sf%theta(mIP) xn(2,i)=sf%zeta(nIP) END DO !m END DO !n WRITE(UNIT_stdOut,'(6X,A)') ' => APPROXIMATION ERROR COMPARED TO INPUT POINTS:' WRITE(UNIT_stdOut,'(6X,A,E11.3)')' max(|X1_fourier-X_input|)=',& MAXVAL(ABS(x1_fbase_in%evalDOF_xn(sf%ntheta*sf%nzeta,xn,0,X1_b)/scale_minor_radius & -RESHAPE(sf%X,(/sf%ntheta*sf%nzeta/)))) WRITE(UNIT_stdOut,'(6X,A,E11.3)')' max(|X2_fourier-Y_input|)', & MAXVAL(ABS(x2_fbase_in%evalDOF_xn(sf%ntheta*sf%nzeta,xn,0,X2_b)/scale_minor_radius & -RESHAPE(sf%Y,(/sf%ntheta*sf%nzeta/)))) WRITE(UNIT_stdOut,'(A)')' ... CONVERT BOUNDARY DONE.' END SUBROUTINE bff_convert_to_modes !=================================================================================================================================== !! deallocate everything !=================================================================================================================================== SUBROUTINE bff_free(sf) ! MODULES IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES CLASS(t_boundaryFromFile), INTENT(INOUT) :: sf !! self !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES !=================================================================================================================================== IF(.NOT.sf%initialized)RETURN SDEALLOCATE(sf%theta) SDEALLOCATE(sf%zeta) SDEALLOCATE(sf%X) SDEALLOCATE(sf%Y) IF(ALLOCATED(sf%nc))THEN CALL sf%nc%free() DEALLOCATE(sf%nc) END IF sf%initialized=.FALSE. END SUBROUTINE bff_free END MODULE MODgvec_boundaryFromFile