Skip to content

Commit

Permalink
Fix cases where we will still assummed the n index was the same as a …
Browse files Browse the repository at this point in the history 
…the toroidal mdoe.
  • Loading branch information
@cianciosa
cianciosa committed Jan 29, 2024
1 parent b107f80 commit 69d2fec
Show file tree
Hide file tree
Showing 10 changed files with 120 additions and 208 deletions.
42 changes: 21 additions & 21 deletions Sources/check3d.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -323,27 +323,27 @@ SUBROUTINE CheckForces(xc, gc)
DO ntype = 1, ndims
DO n = -ntor, ntor
DO m = 0, mpol
icol = icol+1

xc(icol,js) = eps
CALL funct_island
wplus = wtotal

xc(icol,js) = -eps
CALL funct_island
wmins = wtotal
force = -(wplus-wmins)/(2*eps)

CALL ASSERT(l_getwmhd,'L_GETWMHD = FALSE!')

xc(icol,js) = 0

IF (lWrite) THEN
nt1 = ntype
IF (ntype .GT. 3) nt1 = ntype-3
WRITE (3000+iam,100) m, n, ntype, force, force_array(nt1), &
gc1(icol,js)
END IF
icol = icol+1
xc(icol,js) = eps
CALL funct_island
wplus = wtotal

xc(icol,js) = -eps
CALL funct_island
wmins = wtotal
force = -(wplus-wmins)/(2*eps)
CALL ASSERT(l_getwmhd,'L_GETWMHD = FALSE!')

xc(icol,js) = 0

IF (lWrite) THEN
nt1 = ntype
IF (ntype .GT. 3) nt1 = ntype-3
WRITE (3000+iam,100) m, n, ntype, force, force_array(nt1), &
gc1(icol,js)
END IF

END DO
END DO
Expand Down
55 changes: 19 additions & 36 deletions Sources/diagnostics_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,25 +94,6 @@ SUBROUTINE divb(ns_min, ns_max)
f_cos, nsmin, nsmax)
END IF

! DO m = 0, mpol
! DO n = -ntor, ntor
! divbmnsf(m,n,nl:nh) = ohs*(jbsupsmnsh(m,n,nl+1:nh+1) - &
! jbsupsmnsh(m,n,nl:nh)) &
! - m*(jbsupumnch(m,n,nl+1:nh+1) + &
! jbsupumnch(m,n,nl:nh))*0.5 &
! - n*nfp*(jbsupvmnch(m,n,nl+1:nh+1) + &
! jbsupvmnch(m,n,nl:nh))*0.5
! IF (lasym) THEN
! divbmncf(m,n,nl:nh) = ohs*(jbsupsmnch(m,n,nl+1:nh+1) - &
! jbsupsmnch(m,n,nl:nh)) &
! - m*(jbsupumnsh(m,n,nl+1:nh+1) + &
! jbsupumnsh(m,n,nl:nh))*0.5 &
! + n*nfp*(jbsupvmnsh(m,n,nl+1:nh+1) + &
! jbsupvmnsh(m,n,nl:nh))*0.5
! END IF
! END DO
! END DO

tnorm = hs_i*SUM(jbsupumnch(:,:,nsmin:nsmax)**2 &
+ jbsupvmnch(:,:,nsmin:nsmax)**2 &
+ jbsupsmnsh(:,:,nsmin:nsmax)**2)
Expand Down Expand Up @@ -144,7 +125,7 @@ SUBROUTINE divb(ns_min, ns_max)
!> @param[in] ns_max Upper radial index.
!-------------------------------------------------------------------------------
FUNCTION divdb(ns_min, ns_max)
USE island_params, ONLY: ohs=>ohs_i
USE island_params, ONLY: ohs=>ohs_i, fourier_context
USE hessian, ONLY: gather_array

! Declare Arguments
Expand All @@ -156,6 +137,7 @@ FUNCTION divdb(ns_min, ns_max)
INTEGER :: istat
INTEGER :: m
INTEGER :: n
INTEGER :: n_mode
INTEGER :: nl
INTEGER :: nh
REAL (dp) :: tnorm
Expand Down Expand Up @@ -184,21 +166,22 @@ FUNCTION divdb(ns_min, ns_max)
divbmnsf = 0
divbmncf = 0

DO m = 0, mpol
DO n = -ntor, ntor
divbmnsf(m,n,nl:nh) = -m*(djbsupumnch(m,n,nl+1:nh+1) + &
djbsupumnch(m,n,nl:nh))*0.5 &
- n*nfp*(djbsupvmnch(m,n,nl+1:nh+1) + &
djbsupvmnch(m,n,nl:nh))*0.5 &
+ ohs*(djbsupsmnsh(m,n,nl+1:nh+1) - &
djbsupsmnsh(m,n,nl:nh))
DO n = -ntor, ntor
n_mode = fourier_context%tor_modes(n)*nfp
DO m = 0, mpol
divbmnsf(m,n,nl:nh) = -m*(djbsupumnch(m,n,nl+1:nh+1) + &
& djbsupumnch(m,n,nl:nh))*0.5 &
& - n_mode*(djbsupvmnch(m,n,nl+1:nh+1) + &
& djbsupvmnch(m,n,nl:nh))*0.5 &
& + ohs*(djbsupsmnsh(m,n,nl+1:nh+1) - &
& djbsupsmnsh(m,n,nl:nh))
IF (lasym) THEN
divbmncf(m,n,nl:nh) = m*(djbsupumnsh(m,n,nl+1:nh+1) + &
djbsupumnsh(m,n,nl:nh))*0.5 &
+ n*nfp*(djbsupvmnsh(m,n,nl+1:nh+1) + &
djbsupvmnsh(m,n,nl:nh))*0.5 &
+ ohs*(djbsupsmnch(m,n,nl+1:nh+1) - &
djbsupsmnch(m,n,nl:nh))
divbmncf(m,n,nl:nh) = m*(djbsupumnsh(m,n,nl+1:nh+1) + &
& djbsupumnsh(m,n,nl:nh))*0.5 &
& + n_mode*(djbsupvmnsh(m,n,nl+1:nh+1) + &
& djbsupvmnsh(m,n,nl:nh))*0.5 &
& + ohs*(djbsupsmnch(m,n,nl+1:nh+1) - &
& djbsupsmnch(m,n,nl:nh))
END IF
END DO
END DO
Expand Down Expand Up @@ -428,7 +411,6 @@ SUBROUTINE SIESTA_PROFILES (iunit, arr_value, label, fsq_total)
DO n = -ntor, ntor
DO m = 0, mpol
m_array(m,n) = m
n_array(m,n) = n
END DO
END DO

Expand All @@ -437,7 +419,8 @@ SUBROUTINE SIESTA_PROFILES (iunit, arr_value, label, fsq_total)
WRITE (p_format,1000) mnmax
WRITE (iunit,p_format) ' ',' ','MPOL--->', m_array
WRITE (p_format,1000) mnmax
WRITE (iunit,p_format) 'RADIUS','RMS','NTOR--->', n_array
WRITE (iunit,p_format) 'RADIUS','RMS','NTOR--->', &
& fourier_context%tor_modes

WRITE (p_format,1002) mnmax
DO js = 2, ns
Expand Down
60 changes: 33 additions & 27 deletions Sources/fourier.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1084,8 +1084,7 @@ FUNCTION fourier_get_index(this, n)

DO i = 0, SIGN(UBOUND(this%tor_modes, 1), n), SIGN(1, n)
IF (this%tor_modes(i) .eq. n) THEN
n = this%tor_modes(i)

n = i
fourier_get_index = .true.
EXIT
END IF
Expand Down Expand Up @@ -1138,7 +1137,6 @@ FUNCTION test_fourier()

four => fourier_class(pol, tor, theta, zeta, nfp, .false., tor_modes)
test_fourier = .true.
DEALLOCATE(tor_modes)

DO n = -tor, tor
DO m = 0, pol
Expand All @@ -1153,7 +1151,7 @@ FUNCTION test_fourier()
CALL four%tomnsp(testij, result, f_cos)
test_fourier = test_fourier .and. &
check(testmn, result, pol, tor, &
'cosine_parity_test')
& tor_modes, 'cosine_parity_test')

! Test sine parity transform.
IF (m .eq. m0 .and. n .eq. n0) THEN
Expand All @@ -1165,7 +1163,7 @@ FUNCTION test_fourier()
CALL four%tomnsp(testij, result, f_sin)
test_fourier = test_fourier .and. &
check(testmn, result, pol, tor, &
'sine_parity_test')
& tor_modes, 'sine_parity_test')

! Test u derivatives of cosine parity.
testmn = 0
Expand All @@ -1175,7 +1173,7 @@ FUNCTION test_fourier()
testmn(m,n) = -m
test_fourier = test_fourier .and. &
check(testmn, result, pol, tor, &
'du_cosine_parity_test')
& tor_modes, 'du_cosine_parity_test')

! Test v derivatives of cosine parity.
testmn = 0
Expand All @@ -1185,7 +1183,7 @@ FUNCTION test_fourier()
testmn(m,n) = -n*nfp
test_fourier = test_fourier .and. &
check(testmn, result, pol, tor, &
'dv_cosine_parity_test')
& tor_modes, 'dv_cosine_parity_test')

! Test u derivatives of sine parity.
testmn = 0
Expand All @@ -1195,7 +1193,7 @@ FUNCTION test_fourier()
testmn(m,n) = m
test_fourier = test_fourier .and. &
check(testmn, result, pol, tor, &
'du_sine_parity_test')
& tor_modes, 'du_sine_parity_test')

! Test v derivatives of sine parity.
testmn = 0
Expand All @@ -1204,8 +1202,8 @@ FUNCTION test_fourier()
CALL four%tomnsp(testij, result, f_cos)
testmn(m,n) = n*nfp
test_fourier = test_fourier .and. &
check(testmn, result, pol, tor, &
'dv_sine_parity_test')
& check(testmn, result, pol, tor, &
& tor_modes, 'dv_sine_parity_test')
END DO
END DO

Expand All @@ -1215,6 +1213,8 @@ FUNCTION test_fourier()
DEALLOCATE(testmn)
DEALLOCATE(result)

DEALLOCATE(tor_modes)

END FUNCTION

!*******************************************************************************
Expand Down Expand Up @@ -1258,39 +1258,45 @@ FUNCTION check_mn(expected, received, m, n, name)
!-------------------------------------------------------------------------------
!> @brief Check all test values.
!>
!> @param[in] expected Expected value for the test.
!> @param[in] received Recieved value for the test.
!> @param[in] name Name of the test.
!> @param[in] expected Expected value for the test.
!> @param[in] received Recieved value for the test.
!> @param[in] mpol Number of poloidal modes.
!> @param[in] ntor Number of toroidal modes.
!> @param[in] tor_modes Toroidal modes.
!> @param[in] name Name of the test.
!-------------------------------------------------------------------------------
FUNCTION check_all(expected, received, mpol, ntor, name)
FUNCTION check_all(expected, received, mpol, ntor, tor_modes, name)

IMPLICIT NONE

! Declare Arguments
LOGICAL :: check_all
REAL (rprec), DIMENSION(:,:), INTENT(in) :: expected
REAL (rprec), DIMENSION(:,:), INTENT(in) :: received
INTEGER, INTENT(in) :: mpol
INTEGER, INTENT(in) :: ntor
CHARACTER (len=*), INTENT(in) :: name
LOGICAL :: check_all
REAL (rprec), DIMENSION(:,:), INTENT(in) :: expected
REAL (rprec), DIMENSION(:,:), INTENT(in) :: received
INTEGER, INTENT(in) :: mpol
INTEGER, INTENT(in) :: ntor
INTEGER, DIMENSION(-ntor:ntor), INTENT(in) :: tor_modes
CHARACTER (len=*), INTENT(in) :: name

! Local Variables
INTEGER :: m
INTEGER :: n
INTEGER :: moff
INTEGER :: noff
INTEGER :: m
INTEGER :: n
INTEGER :: n_mode
INTEGER :: moff
INTEGER :: noff

! Start of executable code.
moff = LBOUND(expected, 1)
noff = LBOUND(expected, 2) + ntor

check_all = .true.

DO m = 0, mpol
DO n = -ntor, ntor

DO n = -ntor, ntor
DO m = 0, mpol
check_all = check_all .and. check(expected(m + moff, n + noff), &
received(m + moff, n + noff), &
m, n, name)
m, tor_modes(n), name)
END DO
END DO

Expand Down
4 changes: 3 additions & 1 deletion Sources/metrics.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,9 +122,10 @@ END SUBROUTINE init_metric_elements
!>
!> @param[in] mpol_in Number of SIESTA poloidal modes.
!> @param[in] ntor_in Number of SIESTA toroidal modes.
!> @param[in] nfp_in Number of field periods.
!> @param[in] tor_modes Toroidal mode numbers.
!-------------------------------------------------------------------------------
SUBROUTINE set_grid_sizes(mpol_in, ntor_in, tor_modes)
SUBROUTINE set_grid_sizes(mpol_in, ntor_in, nfp_in, tor_modes)
USE island_params
USE shared_data, ONLY: mblk_size, ndims, lasym

Expand All @@ -133,6 +134,7 @@ SUBROUTINE set_grid_sizes(mpol_in, ntor_in, tor_modes)
! Declare Arguments
INTEGER, INTENT(IN) :: mpol_in
INTEGER, INTENT(IN) :: ntor_in
INTEGER, INTENT(IN) :: nfp_in
INTEGER, DIMENSION(-ntor_in:ntor_in), INTENT(in) :: tor_modes

! Start of executable code
Expand Down
15 changes: 10 additions & 5 deletions Sources/perturbation.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -184,6 +184,7 @@ FUNCTION getwmhd(p)
INTEGER :: iprint
INTEGER :: mres0
INTEGER :: nres0
INTEGER :: n
INTEGER :: isign1
INTEGER :: icount
INTEGER :: irscale
Expand Down Expand Up @@ -283,13 +284,15 @@ FUNCTION getwmhd(p)
END IF

! Scan in radius to determine primary resonance.
DO nres0 = -ntor, ntor
DO n = -ntor, ntor
! Avoid 0/0 chip/phip at origin
jstart = 3

nres0 = fourier_context%tor_modes(n)

! Look for multiple resonances
DO WHILE (jstart .lt. ns - 1)

IF (.not.FindResonance(mres0, nres0, rad, jstart)) THEN
EXIT
END IF
Expand All @@ -305,8 +308,9 @@ FUNCTION getwmhd(p)
! vanish at endpoints s=0,1

DO irscale = 1, irad_scale
CALL GetResPert(irscale, isign1, mres0, nres0, rad, &
HelPert0, HelPert0A, chip0, phip0, p_width)
CALL GetResPert(irscale, isign1, mres0, n, rad, &
HelPert0, HelPert0A, chip0, phip0, &
p_width)
xc = 0
! FIXME: Why call this again with the same inputs?
wmhd = getwmhd(xc)
Expand Down Expand Up @@ -343,7 +347,7 @@ FUNCTION getwmhd(p)
CALL FLUSH(unit_out)
END IF

CALL GetResPert(irscale, isign1, mres0, nres0, rad, &
CALL GetResPert(irscale, isign1, mres0, n, rad, &
HelPert0, HelPert0A, chip0, phip0, p_width)

wmhd = getwmhd(xc)
Expand Down Expand Up @@ -503,6 +507,7 @@ SUBROUTINE GetResPert(irscale, isign1, mres0, nres0, rad, &
! Local Variables
INTEGER :: js
INTEGER :: istat
INTEGER :: n
REAL (dp) :: rho
REAL (dp) :: rhores
REAL (dp), DIMENSION(:), ALLOCATABLE :: pert_prof
Expand Down
Loading

0 comments on commit 69d2fec

Please sign in to comment.