Skip to content
Snippets Groups Projects
Commit e71cf68a authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann :seedling:
Browse files

add Cyq0_x0 in the shift in sheared kx connection

parent 2970f894
No related branches found
No related tags found
No related merge requests found
...@@ -84,8 +84,7 @@ CONTAINS ...@@ -84,8 +84,7 @@ CONTAINS
parallel_bc, shift_y, Npol, PB_PHASE parallel_bc, shift_y, Npol, PB_PHASE
READ(lu_in,geometry) READ(lu_in,geometry)
PB_PHASE = .false. PB_PHASE = .false.
IF(geom .EQ. 'z-pinch') shear = 0._xp IF(shear .GT. 0._xp) SHEARED = .TRUE.
IF(shear .NE. 0._xp) SHEARED = .true.
SELECT CASE(parallel_bc) SELECT CASE(parallel_bc)
CASE ('dirichlet') CASE ('dirichlet')
CASE ('periodic') CASE ('periodic')
...@@ -132,7 +131,7 @@ CONTAINS ...@@ -132,7 +131,7 @@ CONTAINS
IF(MODULO(FLOOR(Npol),2) .EQ. 0) ERROR STOP "Npol must be odd for s-alpha" IF(MODULO(FLOOR(Npol),2) .EQ. 0) ERROR STOP "Npol must be odd for s-alpha"
call eval_salpha_geometry call eval_salpha_geometry
Cyq0_x0 = 1._xp Cyq0_x0 = 1._xp
CASE('z-pinch') CASE('z-pinch','Z-pinch','zpinch','Zpinch')
CALL speak('Z-pinch geometry') CALL speak('Z-pinch geometry')
call eval_zpinch_geometry call eval_zpinch_geometry
SHEARED = .FALSE. SHEARED = .FALSE.
...@@ -295,7 +294,6 @@ CONTAINS ...@@ -295,7 +294,6 @@ CONTAINS
dBdx(iz,eo) = -hatB(iz,eo) ! LB = 1 dBdx(iz,eo) = -hatB(iz,eo) ! LB = 1
dBdy(iz,eo) = 0._xp dBdy(iz,eo) = 0._xp
dBdz(iz,eo) = 0._xp ! Gene put a factor hatB or 1/hatR in this dBdz(iz,eo) = 0._xp ! Gene put a factor hatB or 1/hatR in this
ENDDO ENDDO
ENDDO ENDDO
...@@ -372,7 +370,7 @@ CONTAINS ...@@ -372,7 +370,7 @@ CONTAINS
! get the real mode number (iky starts at 1 and is shifted from paral) ! get the real mode number (iky starts at 1 and is shifted from paral)
mn_y = iky-1+local_nky_offset mn_y = iky-1+local_nky_offset
! Formula for the shift due to shear after Npol turns ! Formula for the shift due to shear after Npol turns
shift = 2._xp*PI*shear*kyarray(iky)*Npol shift = 2._xp*PI*shear*kyarray(iky)*Npol*Cyq0_x0
DO ikx = 1,total_nkx DO ikx = 1,total_nkx
! Usual formula for shifting indices using that dkx = 2pi*shear*dky/Nexc ! Usual formula for shifting indices using that dkx = 2pi*shear*dky/Nexc
ikx_zBC_L(iky,ikx) = ikx-mn_y*Nexc ikx_zBC_L(iky,ikx) = ikx-mn_y*Nexc
...@@ -383,14 +381,16 @@ CONTAINS ...@@ -383,14 +381,16 @@ CONTAINS
IF( (kxarray(ikx) - shift) .LT. kx_min ) THEN ! outside of the frequ domain IF( (kxarray(ikx) - shift) .LT. kx_min ) THEN ! outside of the frequ domain
! print*,( ((kxarray(ikx) - shift)-kx_min) .LT. EPSILON(kx_min) ) ! print*,( ((kxarray(ikx) - shift)-kx_min) .LT. EPSILON(kx_min) )
! print*, kxarray(ikx)-kx_min, EPSILON(kx_min) ! print*, kxarray(ikx)-kx_min, EPSILON(kx_min)
! IF( (ikx-mn_y*Nexc) .LT. 1 ) THEN ! IF( (ikx-mn_y*Nexc) .LT. 1 ) THEN
SELECT CASE(parallel_bc) SELECT CASE(parallel_bc)
CASE ('dirichlet')! connected to 0 CASE ('dirichlet','disconnected') ! connected to 0
ikx_zBC_L(iky,ikx) = -99 ikx_zBC_L(iky,ikx) = -99
CASE ('periodic') CASE ('periodic')
ikx_zBC_L(iky,ikx) = ikx ikx_zBC_L(iky,ikx) = ikx
CASE ('cyclic')! reroute it by cycling through modes CASE ('cyclic')! reroute it by cycling through modes
ikx_zBC_L(iky,ikx) = MODULO(ikx_zBC_L(iky,ikx)-1,total_nkx)+1 ikx_zBC_L(iky,ikx) = MODULO(ikx_zBC_L(iky,ikx)-1,total_nkx)+1
CASE DEFAULT
ERROR STOP "The parallel BC is not recognized (dirichlet, periodic, cyclic or disconnected)"
END SELECT END SELECT
ENDIF ENDIF
ENDDO ENDDO
...@@ -408,7 +408,7 @@ CONTAINS ...@@ -408,7 +408,7 @@ CONTAINS
! get the real mode number (iky starts at 1 and is shifted from paral) ! get the real mode number (iky starts at 1 and is shifted from paral)
mn_y = iky-1+local_nky_offset mn_y = iky-1+local_nky_offset
! Formula for the shift due to shear after Npol ! Formula for the shift due to shear after Npol
shift = 2._xp*PI*shear*kyarray(iky)*Npol shift = 2._xp*PI*shear*kyarray(iky)*Npol*Cyq0_x0
DO ikx = 1,total_nkx DO ikx = 1,total_nkx
! Usual formula for shifting indices ! Usual formula for shifting indices
ikx_zBC_R(iky,ikx) = ikx+mn_y*Nexc ikx_zBC_R(iky,ikx) = ikx+mn_y*Nexc
...@@ -418,13 +418,15 @@ CONTAINS ...@@ -418,13 +418,15 @@ CONTAINS
IF( (kxarray(ikx) + shift) .GT. kx_max ) THEN ! outside of the frequ domain IF( (kxarray(ikx) + shift) .GT. kx_max ) THEN ! outside of the frequ domain
! IF( (ikx+mn_y*Nexc) .GT. total_nkx ) THEN ! outside of the frequ domain ! IF( (ikx+mn_y*Nexc) .GT. total_nkx ) THEN ! outside of the frequ domain
SELECT CASE(parallel_bc) SELECT CASE(parallel_bc)
CASE ('dirichlet') ! connected to 0 CASE ('dirichlet','disconnected') ! connected to 0
ikx_zBC_R(iky,ikx) = -99 ikx_zBC_R(iky,ikx) = -99
CASE ('periodic') ! connected to itself as for shearless CASE ('periodic') ! connected to itself as for shearless
ikx_zBC_R(iky,ikx) = ikx ikx_zBC_R(iky,ikx) = ikx
CASE ('cyclic') CASE ('cyclic')
! write(*,*) 'check',ikx,iky, kxarray(ikx) + shift, '>', kx_max ! write(*,*) 'check',ikx,iky, kxarray(ikx) + shift, '>', kx_max
ikx_zBC_R(iky,ikx) = MODULO(ikx_zBC_R(iky,ikx)-1,total_nkx)+1 ikx_zBC_R(iky,ikx) = MODULO(ikx_zBC_R(iky,ikx)-1,total_nkx)+1
CASE DEFAULT
ERROR STOP "The parallel BC is not recognized (dirichlet, periodic, cyclic or disconnected)"
END SELECT END SELECT
ENDIF ENDIF
ENDDO ENDDO
...@@ -452,6 +454,7 @@ CONTAINS ...@@ -452,6 +454,7 @@ CONTAINS
! enddo ! enddo
! print*, pb_phase_R ! print*, pb_phase_R
! print*, shift_y ! print*, shift_y
! print*, kx_min, kx_max
! stop ! stop
!!!!!!! Example of maps ('x' means connected to 0 value, in the table it is -99) !!!!!!! Example of maps ('x' means connected to 0 value, in the table it is -99)
! dirichlet connection map BC of the RIGHT boundary (z=pi*Npol-dz) ! dirichlet connection map BC of the RIGHT boundary (z=pi*Npol-dz)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment