siesta_displacement.f90

!*******************************************************************************
!
!  Note separating the Doxygen comment block here so detailed decription is
!  found in the Module not the file.
!
!*******************************************************************************
      MODULE siesta_displacement
      USE island_params, ONLY: hs_i, ns=>ns_i
      USE stel_kinds
      USE stel_constants        
      USE nscalingtools, ONLY: startglobrow, endglobrow
      USE timer_mod, ONLY: time_update_upperv
 
      IMPLICIT NONE
 
      CONTAINS
 
!*******************************************************************************
!  UTILITY SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE update_upperv
      USE shared_data, ONLY: lasym, gc_sup, xc, col_scale
      USE siesta_state, ONLY: clear_field_perts
      USE quantities, ONLY: gvsupsmncf => fsupsmncf, gvsupumnsf => fsupumnsf,  &
                            gvsupvmnsf => fsupvmnsf, gvsupsmnsf => fsupsmnsf,  &
                            gvsupumncf => fsupumncf, gvsupvmncf => fsupvmncf
      USE fourier, ONLY: f_sin, f_cos
 
      IMPLICIT NONE
 
!  local variables
      INTEGER   :: istat
      REAL (dp) :: ton
      REAL (dp) :: toff
 
!  Start of executable code
      CALL second0(ton)
 
      CALL clear_field_perts
 
!  Use components of gc_sup as scratch array. gvsup*mn*f componets point to
!  fsup*mn*f values which in turn point to parts of the gc_sup array.
!  @ref ScaleDisplacement sets gc_sup to the xc array then applies the column
!  factors so the gc_sup array is treated as the xc here.
      CALL scaledisplacement(gc_sup, xc, col_scale)
 
      CALL initdisplacement(gvsupsmncf, gvsupumnsf, gvsupvmnsf, f_cos)
      IF (lasym) THEN
         CALL initdisplacement(gvsupsmnsf, gvsupumncf, gvsupvmncf, f_sin)
      END IF
 
      CALL second0(toff)
      time_update_upperv = time_update_upperv + (toff - ton)
 
      END SUBROUTINE update_upperv
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE initdisplacement(gvsupsmnf, gvsupumnf, gvsupvmnf, iparity)
      USE fourier, ONLY: f_cos, f_sin, f_none, f_sum, m0, m1, m2
      USE quantities, ONLY: jvsupsijf, jvsupuijf, jvsupvijf, mpol, ntor
      USE shared_data, ONLY: l_push_s, l_push_u, l_push_v, l_pedge,     &
                             l_push_edge, col_scale
      USE v3_utilities, ONLY: assert_eq
      USE island_params, ONLY: fourier_context
      USE utilities, ONLY: set_bndy_fouier_m0, set_bndy_full_origin
 
      IMPLICIT NONE
 
!  Declare Arguments
      REAL (dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(inout) :: gvsupsmnf
      REAL (dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(inout) :: gvsupumnf
      REAL (dp), DIMENSION(0:mpol,-ntor:ntor,ns), INTENT(inout) :: gvsupvmnf
      INTEGER, INTENT(in)                                       :: iparity
 
!  local variables
      INTEGER                                                   :: fours
      INTEGER                                                   :: fouruv
      INTEGER                                                   :: fcomb
      INTEGER                                                   :: nsmin
      INTEGER                                                   :: nsmax
 
!  Start of executable code
      nsmin = max(1, startglobrow - 1)
      nsmax = min(endglobrow + 1, ns)
 
      IF (iparity .eq. f_cos) THEN
         fours = f_cos
         fouruv = f_sin
         fcomb = f_none
      ELSE
         fours = f_sin
         fouruv = f_cos
         fcomb = f_sum
      END IF
 
      CALL assert_eq(1, lbound(gvsupsmnf,3), 'LBOUND WRONG IN InitDisplacement')
      CALL assert_eq(ns, ubound(gvsupsmnf,3),                                  &
                     'UBOUND WRONG IN InitDisplacement')
 
      CALL set_bndy_fouier_m0(gvsupsmnf, gvsupumnf, gvsupvmnf, fours)
 
!  Origin boundary conditions (evolve m = 1 F_u, F_s, m = 0 F_v)
      IF (nsmin .eq. 1) THEN
!  FIXME: Check if this should be here. Seems to converge faster if this is
!         commented out.
         CALL set_bndy_full_origin(gvsupsmnf, gvsupumnf, gvsupvmnf)
 
         IF (.not.l_push_s) THEN
            gvsupsmnf(m1,:,1) = 0
         END IF
         IF (.not.l_push_u) THEN
            gvsupumnf(m1,:,1) = 0
         END IF
         IF (.not.l_push_v) THEN
            gvsupvmnf(m0,:,1) = 0
         END IF
      END IF
 
!  Edge boundary conditions (B_tangential = 0 => vsups = 0; J^s = 0 => F_v = 0)
!  Note that F_v and F_u are both ~J^s at edge, so only one can be varied to
!  prevent a dependent row in the Hessian matrix
      IF (nsmax .eq. ns) THEN
         gvsupsmnf(:,:,ns) = 0
         IF (l_pedge) THEN
            gvsupumnf(:,:,ns) = 0
         END IF
         IF (.not.l_push_edge) THEN
            gvsupumnf(:,:,ns) = 0
            gvsupvmnf(:,:,ns) = 0
         END IF
      END IF
 
!  Calculate contravariant (SUP) velocities in real space on full mesh
      CALL fourier_context%toijsp(gvsupsmnf(:,:,nsmin:nsmax),                  &
                                  jvsupsijf(:,:,nsmin:nsmax), fcomb, fours)
      CALL fourier_context%toijsp(gvsupumnf(:,:,nsmin:nsmax),                  &
                                  jvsupuijf(:,:,nsmin:nsmax), fcomb, fouruv)
      CALL fourier_context%toijsp(gvsupvmnf(:,:,nsmin:nsmax),                  &
                                  jvsupvijf(:,:,nsmin:nsmax), fcomb, fouruv)
 
      END SUBROUTINE
      
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE scaledisplacement(xc_scratch, xc, colscale)
      USE hessian, ONLY: apply_colscale
      USE shared_data, ONLY: mblk_size
 
      IMPLICIT NONE
 
!  Declare Arguments
      REAL (dp), DIMENSION(mblk_size,ns), INTENT(out) :: xc_scratch
      REAL (dp), DIMENSION(mblk_size,ns), INTENT(in)  :: xc
      REAL (dp), DIMENSION(mblk_size,ns), INTENT(in)  :: colscale
 
!  local variables
      INTEGER                                         :: nsmin
      INTEGER                                         :: nsmax
 
!  Start of executable code
      nsmin = max(1, startglobrow - 1)
      nsmax = min(endglobrow + 1, ns)
 
      xc_scratch(:,nsmin:nsmax) = xc(:,nsmin:nsmax)
      CALL apply_colscale(xc_scratch, colscale, nsmin, nsmax)
 
      END SUBROUTINE
 
      END MODULE