fBase_evalDOF_IP_tens Function

private function fBase_evalDOF_IP_tens(sf, deriv, DOFs) result(y_IP)

evaluate all modes at all interpolation points, making use of the tensor product: 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) :: deriv

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

array of all modes (sf%modes)

Return Value real(kind=wp), (sf%mn_IP)


Calls

proc~~fbase_evaldof_ip_tens~~CallsGraph proc~fbase_evaldof_ip_tens t_fBase%fBase_evalDOF_IP_tens __dgemm_nn __dgemm_nn proc~fbase_evaldof_ip_tens->__dgemm_nn proc~fbase_evaldof_xn t_fBase%fBase_evalDOF_xn proc~fbase_evaldof_ip_tens->proc~fbase_evaldof_xn __matvec_n __matvec_n proc~fbase_evaldof_xn->__matvec_n proc~fbase_eval_xn t_fBase%fBase_eval_xn proc~fbase_evaldof_xn->proc~fbase_eval_xn

Called by

proc~~fbase_evaldof_ip_tens~~CalledByGraph proc~fbase_evaldof_ip_tens t_fBase%fBase_evalDOF_IP_tens proc~base_evaldof t_base%base_evalDOF proc~base_evaldof->proc~fbase_evaldof_ip_tens proc~base_evaldof_all t_base%base_evalDOF_all proc~base_evaldof_all->proc~fbase_evaldof_ip_tens proc~fbase_test fBase_test proc~fbase_test->proc~fbase_evaldof_ip_tens proc~fbase_init t_fBase%fBase_init proc~fbase_test->proc~fbase_init proc~get_boozer_sinterp t_sfl_boozer%Get_Boozer_sinterp proc~get_boozer_sinterp->proc~fbase_evaldof_ip_tens proc~fbase_new fBase_new proc~get_boozer_sinterp->proc~fbase_new proc~get_field Get_Field proc~get_field->proc~fbase_evaldof_ip_tens proc~initaverageaxis InitAverageAxis proc~initaverageaxis->proc~fbase_evaldof_ip_tens proc~lambda_solve Lambda_solve proc~lambda_solve->proc~fbase_evaldof_ip_tens proc~transform_angles_sinterp Transform_Angles_sinterp proc~transform_angles_sinterp->proc~fbase_evaldof_ip_tens proc~transform_angles_sinterp->proc~fbase_new proc~base_test Base_test proc~base_test->proc~base_evaldof proc~buildtransform_sfl t_transform_sfl%BuildTransform_SFL proc~buildtransform_sfl->proc~get_boozer_sinterp proc~buildtransform_sfl->proc~transform_angles_sinterp proc~evalaux EvalAux proc~evalaux->proc~base_evaldof proc~evalaux->proc~base_evaldof_all proc~fbase_init->proc~fbase_test proc~get_boozer get_boozer proc~get_boozer->proc~get_boozer_sinterp proc~gvec_to_jorek_prepare gvec_to_jorek_prepare proc~gvec_to_jorek_prepare->proc~get_field proc~init_la_from_solution Init_LA_from_Solution proc~init_la_from_solution->proc~lambda_solve proc~initsolution InitSolution proc~initsolution->proc~initaverageaxis proc~base_new Base_new proc~base_new->proc~base_test proc~base_new->proc~fbase_new proc~evalenergy EvalEnergy proc~evalenergy->proc~evalaux proc~evalforce EvalForce proc~evalforce->proc~evalaux proc~evaltotals EvalTotals proc~evaltotals->proc~evalaux proc~fbase_copy t_fBase%fBase_copy proc~fbase_copy->proc~fbase_init proc~fbase_new->proc~fbase_init proc~initsolutionmhd3d t_functional_mhd3d%InitSolutionMHD3D proc~initsolutionmhd3d->proc~init_la_from_solution proc~initsolutionmhd3d->proc~initsolution proc~initsolutionmhd3d->proc~evalenergy proc~initsolutionmhd3d->proc~evalforce interface~writestate WriteState proc~initsolutionmhd3d->interface~writestate proc~minimizemhd3d_descent MinimizeMHD3D_descent proc~minimizemhd3d_descent->proc~evalaux proc~minimizemhd3d_descent->proc~evalenergy proc~minimizemhd3d_descent->proc~evalforce proc~minimizemhd3d_descent->interface~writestate proc~bff_convert_to_modes t_boundaryFromFile%bff_convert_to_modes proc~bff_convert_to_modes->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~base_new proc~init_base->proc~fbase_new proc~initmhd3d t_functional_mhd3d%InitMHD3D proc~initmhd3d->proc~base_new proc~initmhd3d->proc~bff_convert_to_modes proc~minimizemhd3d t_functional_mhd3d%MinimizeMHD3D proc~minimizemhd3d->proc~minimizemhd3d_descent proc~readstatefilefromascii ReadStateFileFromASCII proc~readstatefilefromascii->proc~base_new proc~sfl_boozer_new sfl_boozer_new proc~sfl_boozer_new->proc~fbase_new proc~transform_sfl_init t_transform_sfl%transform_SFL_init proc~transform_sfl_init->proc~base_new proc~writestatetoascii WriteStateToASCII proc~writestatetoascii->proc~evaltotals program~gvec_post GVEC_POST program~gvec_post->proc~evalenergy program~gvec_post->proc~evalforce interface~readstate ReadState interface~readstate->proc~readstatefilefromascii 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 interface~writestate->proc~writestatetoascii 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 proc~init_gvec_to_jorek->interface~readstate 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_IP_tens(sf,deriv,DOFs) RESULT(y_IP)
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_fBase), INTENT(IN   ) :: sf     !! self
  INTEGER       , INTENT(IN   ) :: deriv  !! =0: base, =2: dthet , =3: dzeta
  REAL(wp)      , INTENT(IN   ) :: DOFs(:)!! array of all modes (sf%modes)
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp)                      :: y_IP(sf%mn_IP)
!-----------------------------------------------------------------------------------------------------------------------------------
! 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:sf%mn_nyq(1),1:2,-sf%mn_max(2):sf%mn_max(2))
!===================================================================================================================================
  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)
!    DO n=-sf%mn_max(2),sf%mn_max(2)
!      DO i=1,sf%mn_nyq(1)
!        Ctmp(i,1,n)=SUM(sf%base1D_IPthet(i,1,:)*Amn(:,n))
!        Ctmp(i,2,n)=SUM(sf%base1D_IPthet(i,2,:)*Amn(:,n))
!      END DO !i
!    END DO !n
!    k=0
!    DO j=1,sf%mn_nyq(2)
!      DO i=1,sf%mn_nyq(1)
!        k=k+1
!        y_IP(k)=SUM(Ctmp(i,1:2,:)*sf%base1D_IPzeta(1:2,:,j))
!      END DO !i
!    END DO !j
    __DGEMM_NN(Ctmp,2*sf%mn_nyq(1),  mTotal,sf%base1D_IPthet,  mTotal,      nTotal,Amn)
    __DGEMM_NN(y_IP,  sf%mn_nyq(1),2*nTotal,            Ctmp,2*nTotal,sf%mn_nyq(2),sf%base1D_IPzeta)
  CASE(DERIV_THET)
    __DGEMM_NN(Ctmp,2*sf%mn_nyq(1),  mTotal,sf%base1D_dthet_IPthet,  mTotal,      nTotal,Amn)
    __DGEMM_NN(y_IP,  sf%mn_nyq(1),2*nTotal,                  Ctmp,2*nTotal,sf%mn_nyq(2),sf%base1D_IPzeta)
  CASE(DERIV_ZETA)
    __DGEMM_NN(Ctmp,2*sf%mn_nyq(1),  mTotal,sf%base1D_IPthet,  mTotal,      nTotal,Amn)
    __DGEMM_NN(y_IP,  sf%mn_nyq(1),2*nTotal,            Ctmp,2*nTotal,sf%mn_nyq(2),sf%base1D_dzeta_IPzeta)
  CASE DEFAULT  !for other derivatives, resort to not precomputed/ explicit computation:
     y_IP = sf%evalDOF_xn(sf%mn_IP,sf%x_IP,deriv,DOFs)
  END SELECT
END FUNCTION fBase_evalDOF_IP_tens