Skip to content

Commit

Permalink
Merge pull request #60 from ORNL-Fusion/master_dev
Browse files Browse the repository at this point in the history 
Master dev
  • Loading branch information
@cianciosa
cianciosa authored Oct 14, 2024
2 parents 08c1f2b + 31c7d15 commit f9a1933
Show file tree
Hide file tree
Showing 4 changed files with 47 additions and 157 deletions.
1 change: 0 additions & 1 deletion Sources/blocktridiagonalsolver_s.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5297,7 +5297,6 @@ SUBROUTINE GetScaleFactors(js,scalevector)
IF (lequi1) THEN
scalevector(ic) = scalevector(ic) + SUM(coltmp)
IF (js.EQ.endglobrow .AND. js.NE.N) THEN
!WRITE (*,*) iam, js, "Bot"
scalevector(ic) = scalevector(ic) + BotScaleFac(ic)
END IF
ELSE
Expand Down
175 changes: 38 additions & 137 deletions Sources/grid_extension.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,26 +50,12 @@ MODULE grid_extension

REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: r1s_vmec, r1ss_vmec
REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: z1s_vmec, z1ss_vmec

REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: r1su_vmec, z1su_vmec
REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: r1sv_vmec, z1sv_vmec

REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: rmnc_v, zmns_v
REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: rmns_v, zmnc_v

REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: r_ext, z_ext
REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: ru_ext, zu_ext
REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: rv_ext, zv_ext

REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: r1_ves, z1_ves
REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: ru_ves, zu_ves
REAL(dp), ALLOCATABLE, DIMENSION(:,:) :: rv_ves, zv_ves

REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: r_full, z_full
REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: ru_full, zu_full
REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) :: rv_full, zv_full

REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: sqrtg_full

REAL(dp) :: s_edge, ming, maxg

Expand Down Expand Up @@ -97,6 +83,7 @@ SUBROUTINE grid_extender(wout_file, itype, istat)
USE siesta_namelist, ONLY: nsin_ext
USE shared_data, ONLY: lverbose
USE island_params, ONLY: fourier_context
USE metrics, ONLY: LoadGrid

IMPLICIT NONE

Expand Down Expand Up @@ -158,9 +145,7 @@ SUBROUTINE grid_extender(wout_file, itype, istat)
END IF

ALLOCATE(ri(ntheta,nzeta), zi(ntheta,nzeta))
ALLOCATE(r1_ves(ntheta,nzeta), z1_ves(ntheta,nzeta), &
ru_ves(ntheta,nzeta), zu_ves(ntheta,nzeta), &
rv_ves(ntheta,nzeta), zv_ves(ntheta,nzeta))
ALLOCATE(r1_ves(ntheta,nzeta), z1_ves(ntheta,nzeta))

! FIXME: The vessel file could have a different periodicity than the equilbirum.
ALLOCATE(vtor_modes(-ntor_v:ntor_v))
Expand All @@ -180,12 +165,6 @@ SUBROUTINE grid_extender(wout_file, itype, istat)
CALL four%toijsp(rmnc_v, r1_ves, f_none, f_cos)
CALL four%toijsp(zmns_v, z1_ves, f_none, f_sin)

CALL four%toijsp(rmnc_v, ru_ves, f_du, f_cos)
CALL four%toijsp(zmns_v, zu_ves, f_du, f_sin)

CALL four%toijsp(rmnc_v, rv_ves, f_dv, f_cos)
CALL four%toijsp(zmns_v, zv_ves, f_dv, f_sin)

DEALLOCATE(four)

! Compute extended "s" coordinate for s>1, assuming Vp(s) = Vp(s=1)*s
Expand All @@ -207,13 +186,11 @@ SUBROUTINE grid_extender(wout_file, itype, istat)
delv = vol1 - (vol0 + p5*vp1*(s_edge*s_edge - 1))
delv = delv/(s_edge-1)**2

ALLOCATE(r_ext (ntheta,nzeta,2:nsv), z_ext (ntheta,nzeta,2:nsv), &
ru_ext(ntheta,nzeta,2:nsv), zu_ext(ntheta,nzeta,2:nsv), &
rv_ext(ntheta,nzeta,2:nsv), zv_ext(ntheta,nzeta,2:nsv), &
ALLOCATE(r_ext(ntheta,nzeta,2:nsv), &
z_ext(ntheta,nzeta,2:nsv), &
stat=istat)

DO j = 2, nsv

s_ext = 1 + (j-1)*hs_i !at j=nsv, s_ext = s_edge

!V(s) = Vp*(s^2-1)/2 + V0
Expand All @@ -226,10 +203,6 @@ SUBROUTINE grid_extender(wout_file, itype, istat)
IF (interp_type.EQ.LIN) THEN
r_ext(:,:,j) = rho1*r1_vmec + rho*r1_ves
z_ext(:,:,j) = rho1*z1_vmec + rho*z1_ves
ru_ext(:,:,j) = rho1*ru_vmec + rho*ru_ves
zu_ext(:,:,j) = rho1*zu_vmec + rho*zu_ves
rv_ext(:,:,j) = rho1*rv_vmec + rho*rv_ves
zv_ext(:,:,j) = rho1*zv_vmec + rho*zv_ves
ELSE
rho2 = rho*rho
IF (interp_type .eq. QUAD) THEN
Expand All @@ -238,14 +211,6 @@ SUBROUTINE grid_extender(wout_file, itype, istat)
+ r1_vmec*(1-rho2) + r1_ves*rho2
z_ext(:,:,j) = p5*z1s_vmec*(1 - rho2 - rho1*rho1) &
+ z1_vmec*(1-rho2) + z1_ves*rho2
ru_ext(:,:,j) = p5*r1su_vmec*(1 - rho2 - rho1*rho1) &
+ ru_vmec*(1-rho2) + ru_ves*rho2
zu_ext(:,:,j) = p5*z1su_vmec*(1 - rho2 - rho1*rho1) &
+ zu_vmec*(1-rho2) + zu_ves*rho2
rv_ext(:,:,j) = p5*r1sv_vmec*(1 - rho2 - rho1*rho1) &
+ rv_vmec*(1-rho2) + rv_ves*rho2
zv_ext(:,:,j) = p5*z1sv_vmec*(1 - rho2 - rho1*rho1) &
+ zv_vmec*(1-rho2) + zv_ves*rho2
ELSE IF (interp_type .eq. CUB) THEN
rho3 = rho*rho2
CALL ASSERT(.FALSE., 'Unknown interpolation <TYPE> specified')
Expand All @@ -254,23 +219,6 @@ SUBROUTINE grid_extender(wout_file, itype, istat)
END DO

nse = nse - 1 !Not to double count s=1 surface
ALLOCATE(r_full(ntheta,nzeta,nse), z_full(ntheta,nzeta,nse), &
ru_full(ntheta,nzeta,nse), zu_full(ntheta,nzeta,nse), &
rv_full(ntheta,nzeta,nse), zv_full(ntheta,nzeta,nse))

r_full (:,:,1:ns) = r1_i
r_full (:,:,ns + 1:nse) = r_ext
z_full (:,:,1:ns) = z1_i
z_full (:,:,ns + 1:nse) = z_ext
ru_full(:,:,1:ns) = ru_i
ru_full(:,:,ns + 1:nse) = ru_ext
zu_full(:,:,1:ns) = zu_i
zu_full(:,:,ns + 1:nse) = zu_ext
rv_full(:,:,1:ns) = rv_i
rv_full(:,:,ns + 1:nse) = rv_ext
zv_full(:,:,1:ns) = zv_i
zv_full(:,:,ns + 1:nse) = zv_ext

nse = MIN(ns + nsin_ext, nse)

IF (lverbose .and. iam .eq. 0) THEN
Expand All @@ -293,21 +241,8 @@ SUBROUTINE grid_extender(wout_file, itype, istat)
DEALLOCATE(zu_i)
DEALLOCATE(rv_i)
DEALLOCATE(zv_i)
ALLOCATE(r1_i(ntheta,nzeta,ns), z1_i(ntheta,nzeta,ns), &
ru_i(ntheta,nzeta,ns), zu_i(ntheta,nzeta,ns), &
rv_i(ntheta,nzeta,ns), zv_i(ntheta,nzeta,ns))
r1_i = r_full(:,:,1:ns)
z1_i = z_full(:,:,1:ns)
ru_i = ru_full(:,:,1:ns)
zu_i = zu_full(:,:,1:ns)
rv_i = rv_full(:,:,1:ns)
zv_i = zv_full(:,:,1:ns)
END IF

DEALLOCATE(r_full, z_full)
DEALLOCATE(ru_full, zu_full)
DEALLOCATE(rv_full, zv_full)

IF (nsp .NE. ns) THEN
nsp1 = nsp + 1
ALLOCATE(tmp(0:mpol,-ntor:ntor,ns), tmp2(0:mpol,-ntor:ntor, ns))
Expand All @@ -318,10 +253,10 @@ SUBROUTINE grid_extender(wout_file, itype, istat)
zmns_i(0:mpol,-ntor:ntor,ns), stat=istat)
rmnc_i(:,:,1:nsp) = tmp(:,:,1:nsp)
zmns_i(:,:,1:nsp) = tmp2(:,:,1:nsp)
CALL fourier_context%tomnsp(r1_i(:,:,nsp1:ns), rmnc_i(:,:,nsp1:ns), &
f_cos)
CALL fourier_context%tomnsp(z1_i(:,:,nsp1:ns), zmns_i(:,:,nsp1:ns), &
f_sin)
CALL fourier_context%tomnsp(r_ext(:,:,2:nsin_ext + 1), &
rmnc_i(:,:,nsp1:ns), f_cos)
CALL fourier_context%tomnsp(z_ext(:,:,2:nsin_ext + 1), &
zmns_i(:,:,nsp1:ns), f_sin)

! PUT PROFILES ONTO EXTENDED MESH
tmp = 0
Expand Down Expand Up @@ -366,16 +301,10 @@ SUBROUTINE grid_extender(wout_file, itype, istat)
zmnc_i(0:mpol,-ntor:ntor,ns), stat=istat)
rmns_i(:,:,1:nsp) = tmp(:,:,1:nsp)
zmnc_i(:,:,1:nsp) = tmp2(:,:,1:nsp)
CALL fourier_context%tomnsp(r1_i(:,:,nsp1:ns), &
CALL fourier_context%tomnsp(r_ext(:,:,1:nsin_ext), &
& rmns_i(:,:,nsp1:ns), f_sin)
CALL fourier_context%tomnsp(z1_i(:,:,nsp1:ns), &
& zmnc_i(:,:,nsp1:ns), f_cos)
DO js = nsp1, ns
rmns_i(:,:,js) = rmns_i(:,:,js) &
/ fourier_context%orthonorm(0:mpol,-ntor:ntor)
zmnc_i(:,:,js) = zmnc_i(:,:,js) &
/ fourier_context%orthonorm(0:mpol,-ntor:ntor)
END DO
& z_ext(:,:,1:nsin_ext), f_cos)

tmp = 0
tmp(:,:,1:nsp)=jcurrumns
Expand Down Expand Up @@ -404,52 +333,25 @@ SUBROUTINE grid_extender(wout_file, itype, istat)
END DO
END DO

! u = theta, v = zeta
ALLOCATE (r12 (ntheta,nzeta), &
& rs12(ntheta,nzeta), &
& zs12(ntheta,nzeta), &
& ru12(ntheta,nzeta), &
& zu12(ntheta,nzeta), &
& sqrtg_full(ntheta,nzeta,nse), &
& stat = istat)
CALL ASSERT(istat.eq.0,' Allocation error in grid_extender')
sqrtg_full(:,:,1) = 0
DO js = 2, nse
js1 = js-1
r12 = p5*(r1_i(:,:,js) + r1_i(:,:,js1))
rs12 = (r1_i(:,:,js) - r1_i(:,:,js1))*ohs
zs12 = (z1_i(:,:,js) - z1_I(:,:,js1))*ohs
ru12 = p5*(ru_i(:,:,js) + ru_i(:,:,js1))
zu12 = p5*(zu_i(:,:,js) + zu_i(:,:,js1))
sqrtg_full(:,:,js) = r12*(ru12*zs12 - rs12*zu12)
END DO

DEALLOCATE (r12, rs12, zs12, ru12, zu12)

ming = MINVAL(sqrtg_full(:,:,ns:nse))
maxg = MAXVAL(sqrtg_full(:,:,ns:nse))
PRINT *,' IAM: ', iam, 'MING: ', ming,' MAXG: ',maxg
CALL ASSERT(ming*maxg.GT.zero,' Jacobian changed sign in grid_extender')
LoadGrid(istat)

ALLOCATE (vp(nse))
CALL SurfAvgLocal(vp,sqrtg_full,1,nse)
CALL SurfAvgLocal(vp,sqrtg,1,nse)
vp(1) = 0
DO js=2, nse
ming = MINVAL(sqrtg_full(:,:,js) )
maxg = MAXVAL(sqrtg_full(:,:,js))
ming = MINVAL(sqrtg(:,:,js) )
maxg = MAXVAL(sqrtg(:,:,js))
WRITE(5000,'(i4,1p,3E14.4)') js, -vp(js), ming, maxg
END DO
CALL FLUSH(5000)
DEALLOCATE (vp,r_ext, z_ext, ru_ext, zu_ext)
DEALLOCATE (ri,zi)
#endif

CALL LoadGrid(istat)

DEALLOCATE(r1s_vmec)
DEALLOCATE(r1su_vmec)
DEALLOCATE(r1sv_vmec)
DEALLOCATE(z1s_vmec)
DEALLOCATE(z1su_vmec)
DEALLOCATE(z1sv_vmec)

END SUBROUTINE

Expand All @@ -473,20 +375,12 @@ SUBROUTINE init_extender

! First derivatives on s=1 surface
ALLOCATE(r1s_vmec(ntheta,nzeta),z1s_vmec(ntheta,nzeta))
ALLOCATE(r1su_vmec(ntheta,nzeta),z1su_vmec(ntheta,nzeta))
ALLOCATE(r1sv_vmec(ntheta,nzeta),z1sv_vmec(ntheta,nzeta))

fac =.5 !COULD SCAN THIS TO MAKE MAX(sqrt(g) < 0) IF NECESSARY
!! fac =.17 !COULD SCAN THIS TO MAKE MAX(sqrt(g) < 0) IF NECESSARY
r1s_vmec=(r1_vmec - r1_i(:,:,ns-1))*ohs*fac
z1s_vmec=(z1_vmec - z1_i(:,:,ns-1))*ohs*fac

r1su_vmec = (ru_vmec - ru_i(:,:,ns-1))*ohs*fac
z1su_vmec = (zu_vmec - zu_i(:,:,ns-1))*ohs*fac

r1sv_vmec = (rv_vmec - rv_i(:,:,ns-1))*ohs*fac
z1sv_vmec = (zv_vmec - zv_i(:,:,ns-1))*ohs*fac

END SUBROUTINE

!-------------------------------------------------------------------------------
Expand All @@ -499,17 +393,18 @@ FUNCTION GetVP1()
IMPLICIT NONE

! Declare Arguments
REAL(dp) :: GetVP1
REAL(dp) :: GetVP1

! Local Variables
INTEGER :: js
INTEGER :: js1
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: r12
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: rs12
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zs12
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: ru12
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zu12
REAL(dp), DIMENSION(:), ALLOCATABLE :: vp
INTEGER :: js
INTEGER :: js1
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: r12
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: rs12
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zs12
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: ru12
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zu12
REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: sqrtg_full
REAL(dp), DIMENSION(:), ALLOCATABLE :: vp

! Start of executable code.
ALLOCATE(vp(ns - 1:ns), sqrtg_full(ntheta,nzeta,ns - 1:ns))
Expand Down Expand Up @@ -618,14 +513,20 @@ FUNCTION get_volume(rho)
DO i = 1, ntheta-1
zu = (zi(i+1,j)-zi(i,j))*odelu
r1 = p5*(ri(i+1,j)+ri(i,j))
get_volume = get_volume + r1**2 * zu
END DO
IF (lasym) THEN
zu = (zi(1,j)-zi(ntheta,j))*odelu
r1 = p5*(ri(1,j)+ri(ntheta,j))
get_volume = get_volume + r1**2*zu
END IF
END DO
END DO

IF (lasym) THEN
DO j = 1, nzeta
DO i = 1, ntheta-1
zu = (zi(i+1,j)-zi(i,j))*odelu
r1 = p5*(ri(i+1,j)+ri(i,j))
get_volume = get_volume + r1**2*zu
END DO
END DO
END IF

get_volume = dnorm*get_volume

END FUNCTION get_volume
Expand Down
18 changes: 9 additions & 9 deletions Sources/metrics.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -582,17 +582,17 @@ SUBROUTINE tolowerh(xsupsij, xsupuij, xsupvij, &
INTEGER, INTENT(IN) :: nsmax

! Start executable code.
xsubsij = gss(:,nsmin:nsmax)*xsupsij &
+ gsu(:,nsmin:nsmax)*xsupuij &
+ gsv(:,nsmin:nsmax)*xsupvij
xsubsij(:,nsmin:nsmax) = gss(:,nsmin:nsmax)*xsupsij(:,nsmin:nsmax) &
+ gsu(:,nsmin:nsmax)*xsupuij(:,nsmin:nsmax) &
+ gsv(:,nsmin:nsmax)*xsupvij(:,nsmin:nsmax)

xsubuij = gsu(:,nsmin:nsmax)*xsupsij &
+ guu(:,nsmin:nsmax)*xsupuij &
+ guv(:,nsmin:nsmax)*xsupvij
xsubuij(:,nsmin:nsmax) = gsu(:,nsmin:nsmax)*xsupsij(:,nsmin:nsmax) &
+ guu(:,nsmin:nsmax)*xsupuij(:,nsmin:nsmax) &
+ guv(:,nsmin:nsmax)*xsupvij(:,nsmin:nsmax)

xsubvij = gsv(:,nsmin:nsmax)*xsupsij &
+ guv(:,nsmin:nsmax)*xsupuij &
+ gvv(:,nsmin:nsmax)*xsupvij
xsubvij(:,nsmin:nsmax) = gsv(:,nsmin:nsmax)*xsupsij(:,nsmin:nsmax) &
+ guv(:,nsmin:nsmax)*xsupuij(:,nsmin:nsmax) &
+ gvv(:,nsmin:nsmax)*xsupvij(:,nsmin:nsmax)

END SUBROUTINE tolowerh

Expand Down
10 changes: 0 additions & 10 deletions Sources/siesta_force.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,16 +118,6 @@ SUBROUTINE update_force
bsupsijh, bsupuijh, bsupvijh, pijh, f_cos)
END IF

! IF (nsmin .le. 5 .and. nsmax .ge. 5) THEN
! DO n = -ntor, ntor
! DO m = 0, mpol
! IF (jpmnch(m,n,5) .gt. 0.0 .or. djpmnch(m,n,6) .gt. 0.0) THEN
! WRITE (*,*) iam, n, m, nsmin, nsmax, jpmnch(m,n,6), djpmnch(m,n,6)
! END IF
! END DO
! END DO
! END IF

! Update thermal energy (pressure based). pijh contains the jacobian term at
! this point.
WPRES: IF (l_getwmhd) THEN
Expand Down

0 comments on commit f9a1933

Please sign in to comment.