hmap_axisNB_init_params Function

public function hmap_axisNB_init_params(ncfile, nvisu) result(sf)

Uses

  • proc~~hmap_axisnb_init_params~~UsesGraph proc~hmap_axisnb_init_params hmap_axisNB_init_params module~modgvec_fbase MODgvec_fBase proc~hmap_axisnb_init_params->module~modgvec_fbase module~modgvec_io_netcdf MODgvec_IO_NETCDF proc~hmap_axisnb_init_params->module~modgvec_io_netcdf module~modgvec_mpi MODgvec_MPI proc~hmap_axisnb_init_params->module~modgvec_mpi module~modgvec_globals MODgvec_Globals module~modgvec_fbase->module~modgvec_globals module~modgvec_io_netcdf->module~modgvec_globals iso_fortran_env iso_fortran_env module~modgvec_globals->iso_fortran_env

initialize the type hmap_axisNB and read "G-frame" from netcdf

transform from full period cartesian coordinates to one-field period "hat" cartesian coordinates, by rotating around the zaxis with the local angle zeta, with the sign given by which direction xyz(zeta=0) rotates to xyz(zeta=2pi/nfp) INVERSION OF: x=cos(zeta)xhat - sgn_rotsin(zeta)yhat, y=cos(zeta)yhat + sgn_rotsin(zeta)xhat, z <=> zhat ==> xhat=cos(zeta)x+sgn_rotsin(zeta)y, yhat=cos(zeta)y-sgn_rotsin(zeta)x, check that all points on full period are the same in the xhat,yhat,zhat coordinates NOTE THIS FUNCTION IS USING sinz=sin(zeta),cosz=cos(zeta) and sgn_rot computed above!!!!

Arguments

Type IntentOptional Attributes Name
character(len=*), intent(in) :: ncfile

netcdf file containing the group "axis" from which to read the G-frame
integer, intent(in) :: nvisu

number of visualization points for G-Frame per field period, -1: no visualization

Return Value type(t_hmap_axisNB)

self


Calls

proc~~hmap_axisnb_init_params~~CallsGraph proc~hmap_axisnb_init_params hmap_axisNB_init_params interface~cross CROSS proc~hmap_axisnb_init_params->interface~cross interface~par_bcast par_Bcast proc~hmap_axisnb_init_params->interface~par_bcast proc~allocate_readin_vars allocate_readin_vars proc~hmap_axisnb_init_params->proc~allocate_readin_vars proc~checkfieldperiodicity CheckFieldPeriodicity proc~hmap_axisnb_init_params->proc~checkfieldperiodicity proc~fbase_initdof t_fBase%fBase_initDOF proc~hmap_axisnb_init_params->proc~fbase_initdof proc~hmap_axisnb_eval_tnb_hat t_hmap_axisNB%hmap_axisNB_eval_TNB_hat proc~hmap_axisnb_init_params->proc~hmap_axisnb_eval_tnb_hat proc~hmap_axisnb_test hmap_axisNB_test proc~hmap_axisnb_init_params->proc~hmap_axisnb_test proc~ncfile_init ncfile_init proc~hmap_axisnb_init_params->proc~ncfile_init proc~par_barrier par_Barrier proc~hmap_axisnb_init_params->proc~par_barrier proc~readnetcdf~2 ReadNETCDF proc~hmap_axisnb_init_params->proc~readnetcdf~2 proc~rodrigues rodrigues proc~hmap_axisnb_init_params->proc~rodrigues proc~visu_axisnb Visu_axisNB proc~hmap_axisnb_init_params->proc~visu_axisnb interface~cross->interface~cross proc~par_bcast_array1d par_Bcast_array1D interface~par_bcast->proc~par_bcast_array1d proc~par_bcast_array1d_int par_Bcast_array1D_int interface~par_bcast->proc~par_bcast_array1d_int proc~par_bcast_array1d_str par_Bcast_array1D_str interface~par_bcast->proc~par_bcast_array1d_str proc~par_bcast_array2d par_Bcast_array2D interface~par_bcast->proc~par_bcast_array2d proc~par_bcast_scalar par_Bcast_scalar interface~par_bcast->proc~par_bcast_scalar proc~par_bcast_scalar_int par_Bcast_scalar_int interface~par_bcast->proc~par_bcast_scalar_int proc~par_bcast_scalar_str par_Bcast_scalar_str interface~par_bcast->proc~par_bcast_scalar_str proc~checkfieldperiodicity->proc~hmap_axisnb_eval_tnb_hat proc~checkfieldperiodicity->proc~rodrigues proc~fbase_projectiptodof_tens t_fBase%fBase_projectIPtoDOF_tens proc~fbase_initdof->proc~fbase_projectiptodof_tens proc~fbase_projectxntodof t_fBase%fBase_projectxntoDOF proc~fbase_initdof->proc~fbase_projectxntodof dgemv dgemv proc~hmap_axisnb_eval_tnb_hat->dgemv proc~fbase_eval t_fBase%fBase_eval proc~hmap_axisnb_eval_tnb_hat->proc~fbase_eval proc~hmap_axisnb_eval t_hmap_axisNB%hmap_axisNB_eval proc~hmap_axisnb_test->proc~hmap_axisnb_eval proc~hmap_axisnb_eval_all t_hmap_axisNB%hmap_axisNB_eval_all proc~hmap_axisnb_test->proc~hmap_axisnb_eval_all proc~hmap_axisnb_eval_dxdq t_hmap_axisNB%hmap_axisNB_eval_dxdq proc~hmap_axisnb_test->proc~hmap_axisnb_eval_dxdq proc~hmap_axisnb_eval_gij t_hmap_axisNB%hmap_axisNB_eval_gij proc~hmap_axisnb_test->proc~hmap_axisnb_eval_gij proc~hmap_axisnb_eval_gij_dq t_hmap_axisNB%hmap_axisNB_eval_gij_dq proc~hmap_axisnb_test->proc~hmap_axisnb_eval_gij_dq proc~hmap_axisnb_eval_jh t_hmap_axisNB%hmap_axisNB_eval_Jh proc~hmap_axisnb_test->proc~hmap_axisnb_eval_jh proc~hmap_axisnb_eval_jh_dq t_hmap_axisNB%hmap_axisNB_eval_Jh_dq proc~hmap_axisnb_test->proc~hmap_axisnb_eval_jh_dq proc~hmap_axisnb_get_ddx_dqij t_hmap_axisNB%hmap_axisNB_get_ddx_dqij proc~hmap_axisnb_test->proc~hmap_axisnb_get_ddx_dqij proc~hmap_axisnb_get_dx_dqi t_hmap_axisNB%hmap_axisNB_get_dx_dqi proc~hmap_axisnb_test->proc~hmap_axisnb_get_dx_dqi proc~hmap_axisnb_init_aux hmap_axisNB_init_aux proc~hmap_axisnb_test->proc~hmap_axisnb_init_aux proc~mpi_check_single_access mpi_check_single_access proc~ncfile_init->proc~mpi_check_single_access proc~ncfile_openfile t_ncfile%ncfile_openfile proc~ncfile_init->proc~ncfile_openfile proc~readnetcdf~2->proc~allocate_readin_vars proc~ncfile_closefile t_ncfile%ncfile_closefile proc~readnetcdf~2->proc~ncfile_closefile proc~ncfile_get_array t_ncfile%ncfile_get_array proc~readnetcdf~2->proc~ncfile_get_array proc~ncfile_get_scalar t_ncfile%ncfile_get_scalar proc~readnetcdf~2->proc~ncfile_get_scalar proc~rodrigues->interface~cross proc~visu_axisnb->proc~hmap_axisnb_eval_tnb_hat proc~writedatatovtk WriteDataToVTK proc~visu_axisnb->proc~writedatatovtk writedatatocsv writedatatocsv proc~visu_axisnb->writedatatocsv proc~fbase_eval_xn t_fBase%fBase_eval_xn proc~fbase_eval->proc~fbase_eval_xn dgemm dgemm proc~fbase_projectiptodof_tens->dgemm proc~fbase_projectxntodof->dgemv proc~fbase_projectxntodof->proc~fbase_eval_xn proc~hmap_axisnb_eval->proc~hmap_axisnb_eval_tnb_hat proc~hmap_axisnb_eval_all_e hmap_axisNB_eval_all_e proc~hmap_axisnb_eval_all->proc~hmap_axisnb_eval_all_e proc~hmap_axisnb_eval_dxdq->proc~hmap_axisnb_eval_tnb_hat proc~hmap_axisnb_eval_gij->proc~hmap_axisnb_eval_tnb_hat proc~hmap_axisnb_eval_gij_dq->proc~hmap_axisnb_eval_tnb_hat proc~hmap_axisnb_eval_jh->interface~cross proc~hmap_axisnb_eval_jh->proc~hmap_axisnb_eval_tnb_hat proc~hmap_axisnb_eval_jh_dq->interface~cross proc~hmap_axisnb_eval_jh_dq->proc~hmap_axisnb_eval_tnb_hat proc~hmap_axisnb_get_ddx_dqij->proc~hmap_axisnb_eval_tnb_hat proc~hmap_axisnb_get_dx_dqi->proc~hmap_axisnb_eval_tnb_hat proc~hmap_axisnb_init_aux->interface~cross proc~hmap_axisnb_init_aux->proc~hmap_axisnb_eval_tnb_hat proc~ncfile_closefile->proc~mpi_check_single_access proc~ncfile_get_array->proc~mpi_check_single_access proc~ncfile_enter_groups t_ncfile%ncfile_enter_groups proc~ncfile_get_array->proc~ncfile_enter_groups proc~ncfile_get_scalar->proc~mpi_check_single_access proc~ncfile_get_scalar->proc~ncfile_enter_groups proc~ncfile_openfile->proc~mpi_check_single_access interface~getfreeunit GETFREEUNIT proc~writedatatovtk->interface~getfreeunit interface~getfreeunit->interface~getfreeunit proc~ncfile_enter_groups->proc~mpi_check_single_access proc~ncfile_enter_groups->proc~ncfile_openfile

Called by

proc~~hmap_axisnb_init_params~~CalledByGraph proc~hmap_axisnb_init_params hmap_axisNB_init_params 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~hmap_axisnb_init->proc~hmap_axisnb_init_params

Source Code

FUNCTION hmap_axisNB_init_params(ncfile,nvisu) RESULT(sf)
! MODULES
  USE MODgvec_fbase      ,ONLY: t_fBase
  USE MODgvec_io_netcdf  ,ONLY: ncfile_init
  USE MODgvec_MPI        ,ONLY: par_BCast,par_barrier
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  CHARACTER(LEN=*),INTENT(IN) :: ncfile  !! netcdf file containing the group "axis" from which to read the G-frame
  INTEGER         ,INTENT(IN) :: nvisu   !! number of visualization points for G-Frame per field period, -1: no visualization
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  TYPE(t_hmap_axisNB)  :: sf !! self
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER :: i
  INTEGER :: error_nfp
  REAL(wp) :: Jh,min_Jh,max_Jh
  REAL(wp),DIMENSION(3) :: X0,T,N,B,Np,Bp
  REAL(wp),ALLOCATABLE :: cosz(:),sinz(:)
!===================================================================================================================================
  CALL par_Barrier(beforeScreenOut='INIT HMAP :: axisNB FRAME OF A CLOSED CURVE ...')

  sf%ncfile=ncfile
  IF(MPIroot)THEN
    ! read axis from netcdf
    CALL ncfile_init(sf%nc,sf%ncfile,"r")
    CALL ReadNETCDF(sf)


    sf%sgn_rot=1
    IF(sf%nfp.GT.2)THEN !no need to check for nfp=1 and nfp=2 is either a positive or negative rotation by pi
      IF(SUM((sf%xyz(:,1+sf%nzeta) -  rodrigues(sf%xyz(:,1),TWOPI/REAL(sf%nfp,wp)))**2).LT.1.0e-12)THEN
         sf%sgn_rot=1
      ELSEIF(SUM((sf%xyz(:,1+sf%nzeta) -  rodrigues(sf%xyz(:,1),-TWOPI/REAL(sf%nfp,wp)))**2).LT.1.0e-12)THEN
         sf%sgn_rot=-1
      ELSE
        CALL abort(__STAMP__, &
        'hmap_axisNB(G-Frame): problem with check of field period rotation: point at zeta=0 does not rotate to point at zeta= +/-2pi/nfp.', &
           TypeInfo="InitializationError")
      END IF
    END IF
    WRITE(UNIT_stdOut,'(4X,A,I2)')'INFO: sign of the rotation from zeta=0 to zeta=2pi/nfp is: ',sf%sgn_rot

    sf%n_max=(sf%nzeta-1)/2 ! maximum mode number on a field period

  END IF !MPIroot
  CALL par_BCast(sf%nzeta,0)
  CALL par_BCast(sf%nfp,0)
  CALL par_BCast(sf%n_max,0)
  CALL par_Bcast(sf%sgn_rot,0)
  IF(.NOT.MPIroot) CALL allocate_readin_vars(sf)
  CALL par_Bcast(sf%zeta,0)
  CALL par_Bcast(sf%xyz,0)
  CALL par_Bcast(sf%Nxyz,0)
  CALL par_Bcast(sf%Bxyz,0)
  !Fourier 1D base on one field period for "hat" coordinates
  sf%fb_hat = t_fBase((/0,sf%n_max/),(/1,sf%nzeta/),sf%nfp,"_sincos_",.FALSE.)
  ALLOCATE(sf%xyz_hat_modes( 3,sf%fb_hat%modes),&
           sf%Nxyz_hat_modes(3,sf%fb_hat%modes),&
           sf%Bxyz_hat_modes(3,sf%fb_hat%modes))
  IF(MPIroot)THEN
    ALLOCATE(cosz(sf%nzeta*sf%nfp),sinz(sf%nzeta*sf%nfp))
    !build cos(zeta),sin(zeta) on full torus
    DO i=1,sf%nfp
      cosz(sf%nzeta*(i-1)+1:sf%nzeta*i)=COS(sf%zeta(1:sf%nzeta)+TWOPI*(REAL(i-1,wp)/REAL(sf%nfp,wp)))
      sinz(sf%nzeta*(i-1)+1:sf%nzeta*i)=SIN(sf%zeta(1:sf%nzeta)+TWOPI*(REAL(i-1,wp)/REAL(sf%nfp,wp)))
    END DO
    sf%xyz_hat_modes =transform_to_hat(sf%xyz, "xyz to xyzhat"  )
    sf%Nxyz_hat_modes=transform_to_hat(sf%Nxyz,"Nxyz to Nxyzhat")
    sf%Bxyz_hat_modes=transform_to_hat(sf%Bxyz,"Bxyz to Bxyzhat")
    DEALLOCATE(cosz,sinz)

    IF(nvisu.GT.0) CALL Visu_axisNB(sf,nvisu*sf%nfp)

    CALL CheckFieldPeriodicity(sf,sf%sgn_rot,error_nfp)
    IF(error_nfp.LT.0) &
       CALL abort(__STAMP__, &
          "hmap_axisNB(G-Frame): check Field Periodicity failed!", &
           TypeInfo="InitializationError")
  END IF !MPIroot

  CALL par_BCast(sf%xyz_hat_modes,0)
  CALL par_BCast(sf%Nxyz_hat_modes,0)
  CALL par_BCast(sf%Bxyz_hat_modes,0)

  ! check right-handedness of T.(N X B ) >0 at X0(zeta) points
  min_Jh= HUGE(1.0_wp)
  max_Jh=-HUGE(1.0_wp)
  DO i=1,97
    CALL sf%eval_TNB(TWOPI*REAL(i-1,wp)/(sf%nfp*97.0_wp),X0,T,N,B,Np,Bp)
    Jh=SUM(T*(CROSS(N,B)))
    min_Jh = MIN(min_Jh,Jh)
    max_Jh = MAX(max_Jh,Jh)
  END DO
  IF(min_Jh.LT.0.0_wp)THEN
    IF(max_Jh.LT.0.0_wp)THEN
      CALL abort(__STAMP__, &
       "hmap_axisNB(G-Frame): Jacobian evaluated at X0: negative everywhere, left-handed coordinates", &
       TypeInfo="InitializationError")
    ELSE
      CALL abort(__STAMP__, &
       "hmap_axisNB(G-Frame): Jacobian evaluated at X0: positive and negative values found!", &
       TypeInfo="InitializationError")
    END IF
  END IF

  sf%initialized=.TRUE.
  CALL par_barrier(afterScreenOut='...DONE')

  IF(.NOT.test_called) CALL hmap_axisNB_test(sf)

  !------------------
  CONTAINS
  !------------------
  !! transform from full period cartesian coordinates to one-field period "hat" cartesian coordinates,
  !! by rotating around the zaxis with the local angle zeta,
  !! with the sign given by which direction xyz(zeta=0) rotates to xyz(zeta=2pi/nfp)
  !! INVERSION OF: x=cos(zeta)*xhat - sgn_rot*sin(zeta)*yhat,
  !!               y=cos(zeta)*yhat + sgn_rot*sin(zeta)*xhat, z <=> zhat
  !! ==>           xhat=cos(zeta)x+sgn_rot*sin(zeta)y,
  !!               yhat=cos(zeta)y-sgn_rot*sin(zeta)x,
  !! check that all points on full period are the same in the xhat,yhat,zhat coordinates
  !! NOTE THIS FUNCTION IS USING sinz=sin(zeta),cosz=cos(zeta) and sgn_rot computed above!!!!
  FUNCTION transform_to_hat(xyz_in,msg)  RESULT(to_hat_modes)
    IMPLICIT NONE

    REAL(wp),INTENT(IN) :: xyz_in(3,sf%nzeta*sf%nfp)
    CHARACTER(*),INTENT(IN) :: msg
    REAL(wp) :: to_hat_modes(3,sf%fb_hat%modes)
    !------------------
    REAL(wp) :: to_hat(3,sf%nzeta*sf%nfp)
    REAL(wp) :: check_xhat(1:3,2:sf%nfp)
    !------------------
    to_hat(1,:)=xyz_in(1,:)*cosz(:)+sf%sgn_rot*xyz_in(2,:)*sinz(:)
    to_hat(2,:)=xyz_in(2,:)*cosz(:)-sf%sgn_rot*xyz_in(1,:)*sinz(:)
    to_hat(3,:)=xyz_in(3,:)
    !check periodicity for remaining nfp
    DO i=2,sf%nfp
      check_xhat(1,i)=MAXVAL(ABS(to_hat(1,1:sf%nzeta)-to_hat(1,sf%nzeta*(i-1)+1:sf%nzeta*i)))
      check_xhat(2,i)=MAXVAL(ABS(to_hat(2,1:sf%nzeta)-to_hat(2,sf%nzeta*(i-1)+1:sf%nzeta*i)))
      check_xhat(3,i)=MAXVAL(ABS(to_hat(3,1:sf%nzeta)-to_hat(3,sf%nzeta*(i-1)+1:sf%nzeta*i)))
    END DO
    IF(ANY(check_xhat.GT.1.0e-12))THEN
      DO i=2,sf%nfp
        WRITE(UNIT_stdout,'(A,I4,A,3E9.2)')'fp=1 vs fp=',i,', check_xyz=',check_xhat(1:3,i)
      END DO
      CALL abort(__STAMP__,&
            "hmap_axisNB(G-Frame): transform from cartesian to hat coordinates"//TRIM(msg)//" yields non-field periodic data!", &
            TypeInfo="InitializationError")
    END IF
    DO i=1,3
      to_hat_modes(i,:)=sf%fb_hat%initDOF(to_hat(i,1:sf%nzeta),thet_zeta_start=(/0.,sf%zeta(1)/))
    END DO
  END FUNCTION transform_to_hat

END FUNCTION hmap_axisNB_init_params