VMEC3D_visu

private subroutine VMEC3D_visu(np_in, minmax, only_planes)

Visualize VMEC flux surface data in planes or 3D, number of radial posisiotns fixed to nFluxVMEC+1, only R,Z,Lambda
Arguments

Type	Intent		Name
integer,	intent(in)	::	np_in(3)	(1) #points in s & theta per element,(2:3) #elements in theta,zeta
real(kind=wp),	intent(in)	::	minmax(3,0:1)	minimum /maximum range in s,theta,zeta [0,1]
logical,	intent(in)	::	only_planes
Calls

Help
Called by

Help
Source Code

SUBROUTINE VMEC3D_visu(np_in,minmax,only_planes)
! MODULES
USE MODgvec_Globals,ONLY:TWOPI,PI,UNIT_stdOut
USE MODgvec_VMEC_Readin
USE MODgvec_VMEC_Vars
USE MODgvec_Output_Vars,ONLY:Projectname
USE MODgvec_Output_vtk,     ONLY: WriteDataToVTK
USE MODgvec_Output_CSV,     ONLY: WriteDataToCSV
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
  INTEGER       , INTENT(IN   ) :: np_in(3)     !! (1) #points in s & theta per element,(2:3) #elements in  theta,zeta
  REAL(wp)      , INTENT(IN   ) :: minmax(3,0:1)  !! minimum /maximum range in s,theta,zeta [0,1]
  LOGICAL       , INTENT(IN   ) :: only_planes  !! true: visualize only planes, false:  full 3D
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER            :: i_s,j_s,i_m,i_n,iMode,nplot(3),mn_IP(2)
  INTEGER,PARAMETER  :: nVal=10
  REAL(wp)           :: coord_visu(     3,nFluxVMEC,np_in(1),np_in(3),np_in(2))
  REAL(wp)           :: var_visu(    nVal,nFluxVMEC,np_in(1),np_in(3),np_in(2))
  REAL(wp)           :: thet(np_in(1),np_in(2)),zeta(np_in(3)),R,Z,LA,sinmn(mn_mode),cosmn(mn_mode)
  REAL(wp)           :: xIP(2),theta_star
  REAL(wp)           :: cosmn_nyq,sinmn_nyq,sqrtg,absB
  CHARACTER(LEN=40)  :: VarNames(nVal)          !! Names of all variables that will be written out
  CHARACTER(LEN=255) :: filename
!===================================================================================================================================
  IF(.NOT.MPIroot) RETURN
  IF(only_planes)THEN
    WRITE(UNIT_stdOut,'(A)') 'Start VMEC visu planes...'
  ELSE
    WRITE(UNIT_stdOut,'(A)') 'Start VMEC visu 3D...'
  END IF
  IF((minmax(1,1)-minmax(1,0)).LE.1e-08)THEN
    WRITE(UNIT_stdOut,'(A,F6.3,A,F6.3)') &
     'WARNING visu3D, nothing to visualize since s-range is <=0, s_min= ',minmax(1,0),', s_max= ',minmax(1,1)
    RETURN
  ELSEIF((minmax(2,1)-minmax(2,0)).LE.1e-08)THEN
    WRITE(UNIT_stdOut,'(A,F6.3,A,F6.3)') &
      'WARNING visu3D, nothing to visualize since theta-range is <=0, theta_min= ',minmax(2,0),', theta_max= ',minmax(2,1)
    RETURN
  ELSEIF((minmax(3,1)-minmax(3,0)).LE.1e-08)THEN
    WRITE(UNIT_stdOut,'(A,F6.3,A,F6.3)') &
      'WARNING visu3D, nothing to visualize since zeta-range is <=0, zeta_min= ',minmax(3,0),', zeta_max= ',minmax(3,1)
    RETURN
  END IF
  VarNames(1)="Phi"
  VarNames(2)="iota"
  VarNames(3)="pressure"
  VarNames(4)="lambda"
  VarNames(5)="sqrtG"
  VarNames(6)="|B|"
  VarNames(7)="s"
  VarNames(8)="theta"
  VarNames(9)="zeta"
  VarNames(10)="Mscale"
  mn_IP(1:2)=np_in(2:3)
  ASSOCIATE(n_s=>nFluxVMEC )
  DO i_m=1,mn_IP(1)
    DO j_s=1,np_in(1)
      thet(j_s,i_m)=TWOPI*(minmax(2,0)+(minmax(2,1)-minmax(2,0)) &
                            *REAL((j_s-1)+(i_m-1)*(np_in(1)-1),wp)/REAL((np_in(1)-1)*mn_IP(1),wp))
    END DO !j_s
  END DO
  DO i_n=1,mn_IP(2)
    zeta(i_n)=TWOPI*(minmax(3,0)+(minmax(3,1)-minmax(3,0))*REAL(i_n-1,wp)/REAL(mn_IP(2)-1,wp))
  END DO

  DO i_s=1,n_s
    var_visu(1,i_s,:,:,:)=Phi_Prof(i_s)
    var_visu(2,i_s,:,:,:)=iotaf(i_s)
    var_visu(3,i_s,:,:,:)=presf(i_s)
    var_visu(7,i_s,:,:,:)=SQRT(Phi_prof(i_s)/Phi_prof(n_s)) !=s
    IF(lasym)THEN
      var_visu(10,i_s,:,:,:) = SUM(xm(:)**(4+1)*(Rmnc(:,i_s)**2+Rmns(:,i_s)**2+Zmnc(:,i_s)**2+Zmns(:,i_s)**2))/&  !pexp=4, qexp=1
                               (SUM(xm(:)**(4  )*(Rmnc(:,i_s)**2+Rmns(:,i_s)**2+Zmnc(:,i_s)**2+Zmns(:,i_s)**2))+1.0e-14)
    ELSE
      var_visu(10,i_s,:,:,:) = SUM(xm(:)**(4+1)*(Rmnc(:,i_s)**2+Zmns(:,i_s)**2))/&  !pexp=4, qexp=1
                               (SUM(xm(:)**(4  )*(Rmnc(:,i_s)**2+Zmns(:,i_s)**2))+1.0e-14)
    END IF
  END DO !i_s
  DO i_n=1,mn_IP(2)
    xIP(2)=zeta(i_n)
    DO i_m=1,mn_IP(1)
      DO j_s=1,np_in(1)
        !xIP(1)=thet(j_s,i_m)
        DO i_s=1,n_s
          !SFL
          XIP(1)=thet(j_s,i_m)
          var_visu(  8,i_s,j_s,i_n,i_m) = xIP(1)
          var_visu(  9,i_s,j_s,i_n,i_m) = xIP(2)

          DO iMode=1,mn_mode
            sinmn(iMode)=SIN(xm(iMode)*xIP(1)-xn(iMode)*xIP(2))
            cosmn(iMode)=COS(xm(iMode)*xIP(1)-xn(iMode)*xIP(2))
          END DO !iMode
          R=0.0_wp
          Z=0.0_wp
          LA=0.0_wp
          DO iMode=1,mn_mode
            R =R +Rmnc(iMode,i_s)*cosmn(iMode)
            Z =Z +Zmns(iMode,i_s)*sinmn(iMode)
            LA=LA+Lmns(iMode,i_s)*sinmn(iMode)
          END DO !iMode
          IF(lasym)THEN
            DO iMode=1,mn_mode
              R =R +Rmns(iMode,i_s)*sinmn(iMode)
              Z =Z +Zmnc(iMode,i_s)*cosmn(iMode)
              LA=LA+Lmnc(iMode,i_s)*cosmn(iMode)
            END DO !iMode
          END IF !lasym
          coord_visu(1,i_s,j_s,i_n,i_m) = R*COS(xIP(2))
          coord_visu(2,i_s,j_s,i_n,i_m) =-R*SIN(xIP(2)) !vmec data(R,phi,Z) was flipped in input to (R,Z,phi)!!
          coord_visu(3,i_s,j_s,i_n,i_m) = Z
          var_visu(  4,i_s,j_s,i_n,i_m) = LA
          sqrtg=0.
          absB =0.
          DO iMode=1,mn_mode_nyq
            cosmn_nyq=COS(xm_nyq(iMode)*xIP(1)-xn_nyq(iMode)*xIP(2))
            sqrtg=sqrtg+gmnc(iMode,i_s)*cosmn_nyq
            absB =absB +bmnc(iMode,i_s)*cosmn_nyq
          END DO !iMode
          IF(lasym)THEN
            DO iMode=1,mn_mode_nyq
              sinmn_nyq=SIN(xm_nyq(iMode)*xIP(1)-xn_nyq(iMode)*xIP(2))
              sqrtg=sqrtg+gmns(iMode,i_s)*sinmn_nyq
              absB =absB +bmns(iMode,i_s)*sinmn_nyq
            END DO !iMode
          END IF !lasym
          var_visu(  5,i_s,j_s,i_n,i_m) = sqrtg*2.0*sqrt(phi_Prof(i_s)/phi_prof(n_s))  !VMEC: s=Phi_norm, but should match GVEC s~sqrt(phi_norm)
          var_visu(  6,i_s,j_s,i_n,i_m) = absB
        END DO !i_s=1,n_s
      END DO !j_s=1,np_in(1)
    END DO !i_n
  END DO !i_m
  !overwrite data on axis with theta average of first flux surface (index 2)
  DO i_n=1,mn_IP(2)
    var_visu(4,1,:,i_n,:) =SUM(var_visu(4,2,:,i_n,:))/REAL(mn_IP(1)*np_in(1))
    var_visu(5,1,:,i_n,:) =SUM(var_visu(5,2,:,i_n,:))/REAL(mn_IP(1)*np_in(1))
    var_visu(6,1,:,i_n,:) =SUM(var_visu(6,2,:,i_n,:))/REAL(mn_IP(1)*np_in(1))
  END DO !i_m
  var_visu(10,1,:,:,:)=1.

  !make grid exactly periodic
  !make theta direction exactly periodic
  IF(ABS((minMax(2,1)-minmax(2,0))-1.0_wp).LT.1.0e-04)THEN !fully periodic
    coord_visu( :,:,np_in(1),:,mn_IP(1))=coord_visu(:,:,1,:,1)
  END IF
  !make zeta direction exactly periodic, only for 3Dvisu
  IF(.NOT.only_planes)THEN
    IF(ABS((minMax(3,1)-minmax(3,0))-1.0_wp).LT.1.0e-04)THEN !fully periodic
      coord_visu( :,:,:,mn_IP(2),:)=coord_visu( :,:,:,1,:)
    END IF
  END IF

  IF(only_planes)THEN
    nplot(1:2)=(/n_s,np_in(1) /)-1
    WRITE(filename,'(A,A)')TRIM(Projectname),"_visu_planes_VMEC.vtu"
    CALL WriteDataToVTK(2,3,nVal,nplot(1:2),mn_IP(1)*mn_IP(2),VarNames, &
                        coord_visu(:,:,:,:,:), &
                          var_visu(:,:,:,:,:),TRIM(filename))
  ELSE
    !3D
    nplot(1:3)=(/n_s,np_in(1),mn_IP(2)/)-1
    WRITE(filename,'(A,A)')TRIM(Projectname),"_visu_3D_VMEC.vtu"
    CALL WriteDataToVTK(3,3,nVal,nplot,mn_IP(1),VarNames, &
                        coord_visu(:,:,:,:,:), &
                          var_visu(:,:,:,:,:),TRIM(filename))
  END IF
  WRITE(filename,'(A,A)')TRIM(Projectname),"_profile_1D_VMEC.csv"
  CALL WriteDataToCSV(VarNames(:) ,var_visu(:,:,1,1,1) ,TRIM(filename)  &
                                  ,append_in=.FALSE.)

  END ASSOCIATE

  WRITE(UNIT_stdOut,'(A)') '... DONE.'

END SUBROUTINE VMEC3D_visu
VMEC3D_visu Subroutine

Contents

Source Code

private subroutine VMEC3D_visu(np_in, minmax, only_planes)

Uses

Arguments

Calls

Called by

Source Code
