Transform_Angles_sinterp Subroutine

private subroutine Transform_Angles_sinterp(AB_base_in, A_in, q_base_in, q_in, q_name, q_base_out, q_out, B_in)

Uses

  • proc~~transform_angles_sinterp~~UsesGraph proc~transform_angles_sinterp Transform_Angles_sinterp module~modgvec_base MODgvec_base proc~transform_angles_sinterp->module~modgvec_base module~modgvec_fbase MODgvec_fBase proc~transform_angles_sinterp->module~modgvec_fbase module~modgvec_globals MODgvec_Globals proc~transform_angles_sinterp->module~modgvec_globals module~modgvec_sgrid MODgvec_sGrid proc~transform_angles_sinterp->module~modgvec_sgrid module~modgvec_base->module~modgvec_fbase module~modgvec_base->module~modgvec_globals module~modgvec_base->module~modgvec_sgrid module~modgvec_sbase MODgvec_sBase module~modgvec_base->module~modgvec_sbase module~modgvec_fbase->module~modgvec_globals iso_fortran_env iso_fortran_env module~modgvec_globals->iso_fortran_env module~modgvec_sgrid->module~modgvec_globals module~modgvec_sbase->module~modgvec_globals module~modgvec_sbase->module~modgvec_sgrid module~sll_m_bsplines sll_m_bsplines module~modgvec_sbase->module~sll_m_bsplines module~sll_m_spline_interpolator_1d sll_m_spline_interpolator_1d module~modgvec_sbase->module~sll_m_spline_interpolator_1d module~sll_m_spline_matrix sll_m_spline_matrix module~modgvec_sbase->module~sll_m_spline_matrix module~sll_m_assert sll_m_assert module~sll_m_bsplines->module~sll_m_assert module~sll_m_bsplines_base sll_m_bsplines_base module~sll_m_bsplines->module~sll_m_bsplines_base module~sll_m_bsplines_non_uniform sll_m_bsplines_non_uniform module~sll_m_bsplines->module~sll_m_bsplines_non_uniform module~sll_m_bsplines_uniform sll_m_bsplines_uniform module~sll_m_bsplines->module~sll_m_bsplines_uniform module~sll_m_errors sll_m_errors module~sll_m_bsplines->module~sll_m_errors module~sll_m_working_precision sll_m_working_precision module~sll_m_bsplines->module~sll_m_working_precision module~sll_m_spline_interpolator_1d->module~sll_m_spline_matrix module~sll_m_spline_interpolator_1d->module~sll_m_assert module~sll_m_boundary_condition_descriptors sll_m_boundary_condition_descriptors module~sll_m_spline_interpolator_1d->module~sll_m_boundary_condition_descriptors module~sll_m_spline_interpolator_1d->module~sll_m_bsplines_base module~sll_m_spline_interpolator_1d->module~sll_m_errors module~sll_m_spline_1d sll_m_spline_1d module~sll_m_spline_interpolator_1d->module~sll_m_spline_1d module~sll_m_spline_interpolator_1d->module~sll_m_working_precision module~sll_m_spline_matrix->module~sll_m_errors module~sll_m_spline_matrix_banded sll_m_spline_matrix_banded module~sll_m_spline_matrix->module~sll_m_spline_matrix_banded module~sll_m_spline_matrix_base sll_m_spline_matrix_base module~sll_m_spline_matrix->module~sll_m_spline_matrix_base module~sll_m_spline_matrix_dense sll_m_spline_matrix_dense module~sll_m_spline_matrix->module~sll_m_spline_matrix_dense module~sll_m_spline_matrix->module~sll_m_working_precision module~sll_m_boundary_condition_descriptors->module~sll_m_working_precision module~sll_m_bsplines_base->module~sll_m_assert module~sll_m_bsplines_base->module~sll_m_working_precision module~sll_m_bsplines_non_uniform->module~sll_m_assert module~sll_m_bsplines_non_uniform->module~sll_m_bsplines_base module~sll_m_bsplines_non_uniform->module~sll_m_working_precision module~sll_m_bsplines_uniform->module~sll_m_assert module~sll_m_bsplines_uniform->module~sll_m_bsplines_base module~sll_m_bsplines_uniform->module~sll_m_errors module~sll_m_bsplines_uniform->module~sll_m_working_precision module~sll_m_errors->iso_fortran_env module~sll_m_spline_1d->module~sll_m_assert module~sll_m_spline_1d->module~sll_m_bsplines_base module~sll_m_spline_1d->module~sll_m_working_precision module~sll_m_spline_matrix_banded->iso_fortran_env module~sll_m_spline_matrix_banded->module~sll_m_assert module~sll_m_spline_matrix_banded->module~sll_m_errors module~sll_m_spline_matrix_banded->module~sll_m_spline_matrix_base module~sll_m_spline_matrix_banded->module~sll_m_working_precision module~sll_m_spline_matrix_base->module~sll_m_working_precision module~sll_m_spline_matrix_dense->iso_fortran_env module~sll_m_spline_matrix_dense->module~sll_m_assert module~sll_m_spline_matrix_dense->module~sll_m_errors module~sll_m_spline_matrix_dense->module~sll_m_spline_matrix_base module~sll_m_spline_matrix_dense->module~sll_m_working_precision

Transform a function from VMEC angles q(s,theta,zeta) to new angles q(s,theta,zeta) by projection onto the modes of the new angles: sigma_mn(theta,zeta) using a given in s Here, new angles are theta=theta+A(theta,zeta), zeta=zeta+B(theta,zeta), with A,B periodic functions and zero average and same base Note that in this routine, the integral is transformed back to (theta,zeta) q_mn = iint_0^2pi q(theta,zeta) sigma_mn(theta,zeta) dtheta dzeta = iint_0^2pi q(theta,zeta) sigma_mn(theta,zeta) [(1+dA/dtheta)(1+dB/dzeta)-(dA/dzetadB/dzeta)] dtheta dzeta

Arguments

Type IntentOptional Attributes Name
class(t_base), intent(in) :: AB_base_in
real(kind=wp), intent(in) :: A_in(1:AB_base_in%s%nBase,1:AB_base_in%f%modes)
class(t_base), intent(in) :: q_base_in
real(kind=wp), intent(in) :: q_in(1:q_base_in%s%nBase,1:q_base_in%f%modes)
character(len=*), intent(in) :: q_name
class(t_base), intent(in) :: q_base_out
real(kind=wp), intent(inout) :: q_out(q_base_out%s%nBase,1:q_base_out%f%modes)
real(kind=wp), intent(in), optional :: B_in(1:AB_base_in%s%nBase,1:AB_base_in%f%modes)

Calls

proc~~transform_angles_sinterp~~CallsGraph proc~transform_angles_sinterp Transform_Angles_sinterp dgemv dgemv proc~transform_angles_sinterp->dgemv interface~progressbar ProgressBar proc~transform_angles_sinterp->interface~progressbar proc~fbase_eval t_fBase%fBase_eval proc~transform_angles_sinterp->proc~fbase_eval proc~fbase_evaldof_ip_tens t_fBase%fBase_evalDOF_IP_tens proc~transform_angles_sinterp->proc~fbase_evaldof_ip_tens proc~sbase_evaldof2d_s t_sBase%sBase_evalDOF2D_s proc~transform_angles_sinterp->proc~sbase_evaldof2d_s proc~to_spline_with_bc to_spline_with_BC proc~transform_angles_sinterp->proc~to_spline_with_bc interface~progressbar->interface~progressbar proc~fbase_eval_xn t_fBase%fBase_eval_xn proc~fbase_eval->proc~fbase_eval_xn dgemm dgemm proc~fbase_evaldof_ip_tens->dgemm proc~fbase_evaldof_xn t_fBase%fBase_evalDOF_xn proc~fbase_evaldof_ip_tens->proc~fbase_evaldof_xn proc~sbase_evaldof2d_s->dgemv proc~sbase_eval t_sBase%sBase_eval proc~sbase_evaldof2d_s->proc~sbase_eval proc~sbase_applybctodof_lgm t_sBase%sBase_applyBCtoDOF_LGM proc~to_spline_with_bc->proc~sbase_applybctodof_lgm proc~sbase_initdof t_sBase%sBase_initDOF proc~to_spline_with_bc->proc~sbase_initdof proc~fbase_evaldof_xn->dgemv proc~fbase_evaldof_xn->proc~fbase_eval_xn proc~solve SOLVE proc~sbase_applybctodof_lgm->proc~solve eval_basis eval_basis proc~sbase_eval->eval_basis eval_basis_and_n_derivs eval_basis_and_n_derivs proc~sbase_eval->eval_basis_and_n_derivs lagrangeinterpolationpolys lagrangeinterpolationpolys proc~sbase_eval->lagrangeinterpolationpolys proc~sgrid_find_elem t_sGrid%sGrid_find_elem proc~sbase_eval->proc~sgrid_find_elem compute_interpolant compute_interpolant proc~sbase_initdof->compute_interpolant dgetrf dgetrf proc~solve->dgetrf dgetrs dgetrs proc~solve->dgetrs

Called by

proc~~transform_angles_sinterp~~CalledByGraph proc~transform_angles_sinterp Transform_Angles_sinterp proc~buildtransform_sfl t_transform_sfl%BuildTransform_SFL proc~buildtransform_sfl->proc~transform_angles_sinterp

Source Code

SUBROUTINE Transform_Angles_sinterp(AB_base_in,A_in,q_base_in,q_in,q_name,q_base_out,q_out,B_in)
! MODULES
USE MODgvec_Globals,ONLY: UNIT_stdOut,Progressbar,testlevel
USE MODgvec_base   ,ONLY: t_base,base_new
USE MODgvec_sGrid  ,ONLY: t_sgrid
USE MODgvec_fbase  ,ONLY: t_fbase, sin_cos_map
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CLASS(t_Base),INTENT(IN) :: AB_base_in                                    !< basis of A and B
  REAL(wp)     ,INTENT(IN) :: A_in(1:AB_base_in%s%nBase,1:AB_base_in%f%modes) !< coefficients of thet*=thet+A(s,theta,zeta)
  CLASS(t_Base),INTENT(IN) :: q_base_in                                     !< basis of function f
  REAL(wp)     ,INTENT(IN) :: q_in(1:q_base_in%s%nBase,1:q_base_in%f%modes) !< coefficients of f
  CHARACTER(LEN=*),INTENT(IN):: q_name
  CLASS(t_base),INTENT(IN) :: q_base_out                                    !< new fourier basis of function q in new angles, defined mn_max,mn_nyq
  REAL(wp)    ,INTENT(IN),OPTIONAL :: B_in(1:AB_base_in%s%nBase,1:AB_base_in%f%modes) !< coefficients of zeta*=zeta+B(s,theta,zeta)
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  REAL(wp) ,INTENT(INOUT) :: q_out(q_base_out%s%nBase,1:q_base_out%f%modes)          !< coefficients of q in new angles
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER               :: nBase,i,is,i_mn,mn_IP,mn_max(2),mn_nyq(2)
  REAL(wp)              :: spos,dthet_dzeta
  REAL(wp)              :: check(1:7,q_base_out%s%nBase)
  LOGICAL               :: docheck
  LOGICAL               :: Bpresent
  REAL(wp)              :: A_s(1:AB_base_in%f%modes)
  REAL(wp)              :: B_s(1:AB_base_in%f%modes)
  REAL(wp)              :: q_in_s(1:q_base_in%f%modes)
  REAL(wp), ALLOCATABLE :: q_IP(:),q_m(:,:)   ! q evaluated at spos and all integration points
  REAL(wp), ALLOCATABLE :: f_IP(:)       ! =q*(1+dlambda/dtheta) evaluated at integration points
  REAL(wp), ALLOCATABLE :: modes_IP(:,:) ! mn modes of q evaluated at theta*,zeta* for all integration points
  TYPE(t_fBase)                     :: q_fbase_nyq
  TYPE(t_fBase)                     :: AB_fbase_nyq
  REAL(wp),DIMENSION(:),ALLOCATABLE :: A_IP,dAdthet_IP,B_IP,dBdthet_IP,dBdzeta_IP,dAdzeta_IP
!===================================================================================================================================
  docheck=(testlevel.GT.0)
  Bpresent=PRESENT(B_in)
  mn_max(1:2) =q_base_out%f%mn_max
  mn_nyq(1:2) =q_base_out%f%mn_nyq

  SWRITE(UNIT_StdOut,'(A,I4,3(A,2I6),A,L)')'TRANSFORM '//TRIM(q_name)//' TO NEW ANGLE COORDINATES, nfp=',q_base_in%f%nfp, &
                              ', mn_max_in=',q_base_in%f%mn_max,', mn_max_out=',mn_max,', mn_int=',mn_nyq, ', B_in= ',Bpresent

  __PERFON('transform_angles')
  __PERFON('init')
  !initialize


  !total number of integration points
  mn_IP = q_base_out%f%mn_IP
  nBase = q_base_out%s%nBase

  SWRITE(UNIT_StdOut,*)'        ...Init q_out Base Done'


  !same base for X1, but with new mn_nyq (for pre-evaluation of basis functions)
  q_fbase_nyq =  t_fBase(q_base_in%f%mn_max,  mn_nyq, &
                                q_base_in%f%nfp, &
                    sin_cos_map(q_base_in%f%sin_cos), &
                                q_base_in%f%exclude_mn_zero)
  SWRITE(UNIT_StdOut,*)'        ...Init q_nyq Base Done'
  !same base for lambda, but with new mn_nyq (for pre-evaluation of basis functions)
  AB_fbase_nyq =  t_fBase(AB_base_in%f%mn_max,  mn_nyq, &
                                AB_base_in%f%nfp, &
                    sin_cos_map(AB_base_in%f%sin_cos), &
                                AB_base_in%f%exclude_mn_zero)

  SWRITE(UNIT_StdOut,*)'        ...Init AB_nyq Base Done'

  IF(.NOT.Bpresent) THEN
    ALLOCATE(A_IP(1:mn_IP),dAdthet_IP(1:mn_IP))
  ELSE ! Bpresent
    ALLOCATE(A_IP(1:mn_IP),dAdthet_IP(1:mn_IP),dAdzeta_IP(1:mn_IP),B_IP(1:mn_IP),dBdthet_IP(1:mn_IP),dBdzeta_IP(1:mn_IP))
  END IF !.NOT.Bpresent

  ALLOCATE(f_IP(1:mn_IP),q_IP(1:mn_IP),modes_IP(1:q_base_out%f%modes,1:mn_IP))

  ALLOCATE(q_m(1:q_base_out%f%modes,nBase))
  __PERFOFF('init')

  dthet_dzeta  =q_base_out%f%d_thet*q_base_out%f%d_zeta !integration weights

  CALL ProgressBar(0,nBase)!init
  DO is=1,nBase
    __PERFON('eval_data_s')
    spos=MIN(MAX(1.0e-4_wp,q_base_out%s%s_IP(is)),1.0_wp-1.0e-12_wp) !interpolation points for q_in
    !evaluate q_in at spos
    q_in_s(:)= q_base_in%s%evalDOF2D_s(spos,q_base_in%f%modes,   0,q_in(:,:))
    !evaluate A at spos
    A_s(:)= AB_base_in%s%evalDOF2D_s(spos,AB_base_in%f%modes,   0,A_in(:,:))
    IF(Bpresent)THEN
      B_s(:)     = AB_base_in%s%evalDOF2D_s(spos,AB_base_in%f%modes,   0,B_in(:,:))
    END IF
    __PERFOFF('eval_data_s')
    __PERFON('eval_data_f')
    !evaluate lambda at spos
    ! TEST EXACT CASE: A_s=0.

    IF(.NOT.Bpresent)THEN
      q_IP       =  q_fbase_nyq%evalDOF_IP(         0, q_in_s(  :))
      A_IP       = AB_fbase_nyq%evalDOF_IP(         0, A_s(  :))
      dAdthet_IP = AB_fbase_nyq%evalDOF_IP(DERIV_THET, A_s(  :))
!$OMP PARALLEL DO        &
!$OMP   SCHEDULE(STATIC) DEFAULT(NONE)    &
!$OMP   PRIVATE(i_mn)        &
!$OMP   SHARED(mn_IP,q_base_out,q_IP,A_IP,dAdthet_IP,f_IP,modes_IP)
      DO i_mn=1,mn_IP
        f_IP(i_mn) = q_IP(i_mn)*(1.0_wp + dAdthet_IP(i_mn))
        !evaluate (theta*,zeta*) modes of q_in at (theta,zeta)
        modes_IP(:,i_mn)= q_base_out%f%eval(0,(/q_base_out%f%x_IP(1,i_mn)+A_IP(i_mn),q_base_out%f%x_IP(2,i_mn)/))
      END DO !i_mn=1,mn_IP
!$OMP END PARALLEL DO

    ELSE !Bpresent
      q_IP       =  q_fbase_nyq%evalDOF_IP(         0, q_in_s(  :))
      A_IP       = AB_fbase_nyq%evalDOF_IP(         0, A_s(  :))
      dAdthet_IP = AB_fbase_nyq%evalDOF_IP(DERIV_THET, A_s(  :))
      dAdzeta_IP = AB_fbase_nyq%evalDOF_IP(DERIV_ZETA, A_s(  :))
      B_IP       = AB_fbase_nyq%evalDOF_IP(         0, B_s(  :))
      dBdthet_IP = AB_fbase_nyq%evalDOF_IP(DERIV_THET, B_s(  :))
      dBdzeta_IP = AB_fbase_nyq%evalDOF_IP(DERIV_ZETA, B_s(  :))

!$OMP PARALLEL DO        &
!$OMP   SCHEDULE(STATIC) DEFAULT(NONE)    &
!$OMP   PRIVATE(i_mn)        &
!$OMP   SHARED(mn_IP,q_base_out,q_IP,A_IP,dAdthet_IP,B_IP,dBdthet_IP,dBdzeta_IP,dAdzeta_IP,f_IP,modes_IP)
      DO i_mn=1,mn_IP
        f_IP(i_mn) = q_IP(i_mn)*((1.0_wp + dAdthet_IP(i_mn))*(1.0_wp + dBdzeta_IP(i_mn))-dAdzeta_IP(i_mn)*dBdthet_IP(i_mn))
        !evaluate (theta*,zeta*) modes of q_in at (theta,zeta)
        modes_IP(:,i_mn)= q_base_out%f%eval(0,(/q_base_out%f%x_IP(1,i_mn)+A_IP(i_mn),q_base_out%f%x_IP(2,i_mn)+B_IP(i_mn)/))
      END DO !i_mn=1,mn_IP
!$OMP END PARALLEL DO
    END IF !Bpresent
    __PERFON('project')
    __MATVEC_N(q_m(:,is),modes_IP(:,:),f_IP(:))

    __PERFOFF('project')
    __PERFOFF('eval_data_f')

    !projection: integrate (sum over mn_IP), includes normalization of base!
    !q_m(:,is)=(dthet_dzeta*q_base_out%f%snorm_base(:))*(MATMUL(modes_IP(:,1:mn_IP),f_IP(1:mn_IP)))

    q_m(:,is)=dthet_dzeta*q_base_out%f%snorm_base(:)*q_m(:,is)

    !CHECK at all IP points
    IF(doCheck) THEN
      __PERFON('check')
!      __MATVEC_N(f_IP,q_base_out%f%base_IP,q_m(:,is)) !other points
!      check(6)=MIN(check(6),MINVAL(f_IP))
!      check(7)=MAX(check(7),MAXVAL(f_IP))
      check(6,is)=MINVAL(q_IP)
      check(7,is)=MAXVAL(q_IP)
      !f_IP=MATMUL(q_m(:,is),modes_IP(:,:))
      __MATVEC_T(f_IP,modes_IP,q_m(:,is)) !same points
      check(4,is)=MINVAL(f_IP)
      check(5,is)=MAXVAL(f_IP)


      !f_IP = ABS(f_IP - q_IP)/SUM(ABS(q_IP))*REAL(mn_IP,wp)
      f_IP=ABS(f_ip- q_IP)
      check(1,is)=MINVAL(f_IP)
      check(2,is)=MAXVAL(f_IP)
      check(3,is)=SUM(f_IP)/REAL(mn_IP,wp)
      !WRITE(*,*)'     ------  spos= ',spos
      !WRITE(*,*)'check |q_in-q_out|/(surfavg|q_in|) (min/max/avg)',MINVAL(f_IP),MAXVAL(f_IP),SUM(f_IP)/REAL(mn_IP,wp)
      !WRITE(*,*)'max,min,sum q_out |modes|',MAXVAL(ABS(q_out(is,:))),MINVAL(ABS(q_out(is,:))),SUM(ABS(q_out(is,:)))
      __PERFOFF('check')
    END IF !doCheck

    CALL ProgressBar(is,nBase)
  END DO !is

  CALL to_spline_with_BC(q_base_out,q_m,q_out)

  !finalize
  DEALLOCATE( A_IP, dAdthet_IP)
  IF(Bpresent) DEALLOCATE(dAdzeta_IP, B_IP,dBdthet_IP, dBdzeta_IP )

  DEALLOCATE(modes_IP,q_IP,f_IP,q_m)

  IF(doCheck) THEN
    DO i=4,1,-1
      is=MAX(1,nBase/i)
      WRITE(UNIT_StdOut,'(A,E11.4)')'at rho=',q_base_out%s%s_IP(is)
      WRITE(UNIT_StdOut,'(A,2E21.11)') '   MIN/MAX of input  '//TRIM(q_name)//':', check(6:7,is)
      WRITE(UNIT_StdOut,'(A,2E21.11)') '   MIN/MAX of output '//TRIM(q_name)//':', check(4:5,is)
      WRITE(UNIT_StdOut,'(2A,3E11.4)')'    ERROR of '//TRIM(q_name)//':', &
                                      ' |q_in-q_out| (min/max/avg)',check(1:2,is),check(3,is)
    END DO
  END IF

  SWRITE(UNIT_StdOut,'(A)') '...DONE.'
  __PERFOFF('transform_angles')
END SUBROUTINE Transform_Angles_sinterp