fBase_evalDOF_xn_tens Function

private function fBase_evalDOF_xn_tens(sf, nthet, nzeta, thet, zeta, deriv, DOFs) result(y)

evaluate all modes on a tensor-produc grid (t_i,z_j), making use of the tensor product in the fourier series: y_ij = DOFs_mn * SIN(mt_i - nz_j ) => SIN(mt_i) DOFs_mn COS(nz_j) -COS(mt_i) DOFs_mn SIN(nz_j) y_ij = DOFs_mn * COS(mt_i - nz_j ) => COS(mt_i) DOFs_mn COS(nz_j) +SIN(mt_i) DOFs_mn SIN(nz_j) => a1_im DOFs_mn b1_nj + a2_im DOFs_mn b2_nj can be written as 2 SPECIAL MATMAT operations: c(i,1,n)=a1(i,m) DOFs(m,n) , c(i,2,n) = a2(i,m) DOFs(m,n) => c(i,d,n) = DOT_PROD(a(i,d,1:mmax),DOFs(1:mmax,n)) y(i,j) = c(i,1,n) b1(n,j) + c(i,2,n) b2(n,j) = DOT_PROD(c(i,1:2,1:nmax),b(1:2,1:nmax,j) the 1D ordering in y does not neead a reshape, y(i,j) => y(1:mn_IP), 1D array data can be kept, as it is passed (with its start adress) to DGEMM.

Type Bound

t_fBase

Arguments

Type IntentOptional Attributes Name
class(t_fBase), intent(in) :: sf

self
integer, intent(in) :: nthet

number of points in theta
integer, intent(in) :: nzeta

number of points in zeta
real(kind=wp), intent(in) :: thet(1:nthet)

theta positions
real(kind=wp), intent(in) :: zeta(1:nzeta)

zeta positions
integer, intent(in) :: deriv

=0: base, =2: dthet , =3: dzeta
real(kind=wp), intent(in) :: DOFs(:)

array of all modes

Return Value real(kind=wp), (1:nthet*nzeta)

DOFS evaluated on tensor-product grid,


Calls

proc~~fbase_evaldof_xn_tens~~CallsGraph proc~fbase_evaldof_xn_tens t_fBase%fBase_evalDOF_xn_tens __dgemm_nn __dgemm_nn proc~fbase_evaldof_xn_tens->__dgemm_nn proc~fbase_eval1d_thet fBase_eval1d_thet proc~fbase_evaldof_xn_tens->proc~fbase_eval1d_thet proc~fbase_eval1d_zeta fBase_eval1d_zeta proc~fbase_evaldof_xn_tens->proc~fbase_eval1d_zeta

Called by

proc~~fbase_evaldof_xn_tens~~CalledByGraph proc~fbase_evaldof_xn_tens t_fBase%fBase_evalDOF_xn_tens proc~evaluate_base_tens evaluate_base_tens proc~evaluate_base_tens->proc~fbase_evaldof_xn_tens proc~evaluate_base_tens_all evaluate_base_tens_all proc~evaluate_base_tens_all->proc~fbase_evaldof_xn_tens proc~fbase_test fBase_test proc~fbase_test->proc~fbase_evaldof_xn_tens proc~fbase_init t_fBase%fBase_init proc~fbase_test->proc~fbase_init proc~fbase_init->proc~fbase_test proc~fbase_copy t_fBase%fBase_copy proc~fbase_copy->proc~fbase_init proc~fbase_new fBase_new proc~fbase_new->proc~fbase_init proc~base_new Base_new proc~base_new->proc~fbase_new proc~bff_convert_to_modes t_boundaryFromFile%bff_convert_to_modes proc~bff_convert_to_modes->proc~fbase_new proc~get_boozer_sinterp t_sfl_boozer%Get_Boozer_sinterp proc~get_boozer_sinterp->proc~fbase_new proc~hmap_axisnb_init_params hmap_axisNB_init_params proc~hmap_axisnb_init_params->proc~fbase_new proc~init_base Init_Base proc~init_base->proc~fbase_new proc~init_base->proc~base_new proc~sfl_boozer_new sfl_boozer_new proc~sfl_boozer_new->proc~fbase_new proc~transform_angles_sinterp Transform_Angles_sinterp proc~transform_angles_sinterp->proc~fbase_new interface~t_hmap_axisnb t_hmap_axisNB interface~t_hmap_axisnb->proc~hmap_axisnb_init_params proc~hmap_axisnb_init hmap_axisNB_init interface~t_hmap_axisnb->proc~hmap_axisnb_init proc~buildtransform_sfl t_transform_sfl%BuildTransform_SFL proc~buildtransform_sfl->proc~get_boozer_sinterp proc~buildtransform_sfl->proc~transform_angles_sinterp proc~get_boozer get_boozer proc~get_boozer->proc~get_boozer_sinterp proc~hmap_axisnb_init->proc~hmap_axisnb_init_params proc~init_gvec_to_jorek init_gvec_to_jorek proc~init_gvec_to_jorek->proc~init_base interface~readstate ReadState proc~init_gvec_to_jorek->interface~readstate proc~initmhd3d t_functional_mhd3d%InitMHD3D proc~initmhd3d->proc~base_new proc~initmhd3d->proc~bff_convert_to_modes proc~readstatefilefromascii ReadStateFileFromASCII proc~readstatefilefromascii->proc~base_new proc~transform_sfl_init t_transform_sfl%transform_SFL_init proc~transform_sfl_init->proc~base_new interface~readstate->proc~readstatefilefromascii proc~transform_sfl_new transform_sfl_new proc~transform_sfl_new->proc~transform_sfl_init proc~restartfromstate RestartFromState proc~restartfromstate->interface~readstate

Source Code

FUNCTION fBase_evalDOF_xn_tens(sf,nthet,nzeta,thet,zeta,deriv,DOFs) RESULT(y)
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_fBase), INTENT(IN   ) :: sf     !! self
  INTEGER       , INTENT(IN   ) :: nthet  !! number of points in theta
  INTEGER       , INTENT(IN   ) :: nzeta  !! number of points in zeta
  REAL(wp)      , INTENT(IN   ) :: thet(1:nthet) !! theta positions
  REAL(wp)      , INTENT(IN   ) :: zeta(1:nzeta) !! zeta positions
  INTEGER       , INTENT(IN   ) :: deriv  !! =0: base, =2: dthet , =3: dzeta
  REAL(wp)      , INTENT(IN   ) :: DOFs(:)  !! array of all modes
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp)                      :: y(1:nthet*nzeta)   !! DOFS evaluated on tensor-product grid,
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER                       :: iMode,offset,mTotal,nTotal
  REAL(wp)                      :: Amn(1:sf%mTotal1D,-sf%mn_max(2):sf%mn_max(2))
  REAL(wp)                      :: Ctmp(1:nthet,1:2,-sf%mn_max(2):sf%mn_max(2))
  REAL(wp)                      :: base1D_thet(1:nthet,1:2,1:sf%mTotal1D)
  REAL(wp)                      :: base1D_zeta(1:2,-sf%mn_max(2):sf%mn_max(2),1:nzeta)
!===================================================================================================================================
  IF(SIZE(DOFs,1).NE.sf%modes) CALL abort(__STAMP__, &
         'nDOF not correct when calling fBase_evalDOF_IP_tens' )

  offset=sf%mTotal1D-(sf%mn_max(1)+1) !=0 if sin or cos, =sf%mn_max(1)+1 if sin+cos
  !initialize non existing modes to zero
  Amn(1,-sf%mn_max(2):0)=0.0_wp
  IF(offset.GT.0) Amn(offset+1,-sf%mn_max(2):0)=0.0_wp

  !copy DOFs to  (0:m_max , -n_max:n_max) matrix, careful: Xmn(2,:) has nfp factor!
  DO iMode=sf%sin_range(1)+1,sf%sin_range(2)
    Amn(1+sf%Xmn(1,iMode),sf%Xmn(2,iMode)/sf%nfp)=DOFs(iMode)
  END DO
  DO iMode=sf%cos_range(1)+1,sf%cos_range(2)
    Amn(offset+1+sf%Xmn(1,iMode),sf%Xmn(2,iMode)/sf%nfp)=DOFs(iMode)
  END DO

  mTotal=  sf%mTotal1D
  nTotal=2*sf%mn_max(2)+1 !-n_max:n_nax

  SELECT CASE(deriv)
  CASE(0)
    base1d_thet=fBase_eval1d_thet(sf,0,nthet,thet)
    base1d_zeta=fBase_eval1d_zeta(sf,0,nzeta,zeta)
  CASE(DERIV_THET)
    base1d_thet=fBase_eval1d_thet(sf,1,nthet,thet)
    base1d_zeta=fBase_eval1d_zeta(sf,0,nzeta,zeta)
  CASE(DERIV_ZETA)
    base1d_thet=fBase_eval1d_thet(sf,0,nthet,thet)
    base1d_zeta=fBase_eval1d_zeta(sf,1,nzeta,zeta)
  CASE(DERIV_THET_THET)
    base1d_thet=fBase_eval1d_thet(sf,2,nthet,thet)
    base1d_zeta=fBase_eval1d_zeta(sf,0,nzeta,zeta)
  CASE(DERIV_THET_ZETA)
    base1d_thet=fBase_eval1d_thet(sf,1,nthet,thet)
    base1d_zeta=fBase_eval1d_zeta(sf,1,nzeta,zeta)
  CASE(DERIV_ZETA_ZETA)
    base1d_thet=fBase_eval1d_thet(sf,0,nthet,thet)
    base1d_zeta=fBase_eval1d_zeta(sf,2,nzeta,zeta)
  CASE DEFAULT  !for other derivatives, resort to not precomputed/ explicit computation:
    CALL abort(__STAMP__, &
         "fbase_evalDOF_xn_tens: derivative must be 0,DERIV_THET,_ZETA,_THET_THET,_THET_ZETA,_ZETA_ZETA!")
  END SELECT
  __DGEMM_NN(Ctmp,2*nthet,  mTotal,base1D_thet,  mTotal, nTotal,Amn)
  __DGEMM_NN(y   ,  nthet,2*nTotal,       Ctmp,2*nTotal, nzeta ,base1D_zeta)
END FUNCTION fBase_evalDOF_xn_tens