Skip to content

Commit

Permalink
B(0) fix, Plot cleanup, minor code cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
@jcwright77
jcwright77 committed Sep 25, 2024
1 parent 885b800 commit 641606b
Show file tree
Hide file tree
Showing 7 changed files with 180 additions and 126 deletions.
80 changes: 53 additions & 27 deletions src/CQL3D_SETUP/ceez.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -239,7 +239,7 @@ subroutine curv1 (n,x,y,slp1,slpn,islpsw,yp,temp,
yp(1) = (dx1-slpp1)/diag1
temp(1) = sdiag1/diag1
if (n .eq. 2) go to 6
do 5 i = 2,nm1
do i = 2,nm1
delx2 = x(i+1)-x(i)
if (delx2 .le. 0.) go to 9
dx2 = (y(i+1)-y(i))/delx2
Expand All @@ -249,15 +249,17 @@ subroutine curv1 (n,x,y,slp1,slpn,islpsw,yp,temp,
temp(i) = sdiag2/diag
dx1 = dx2
diag1 = diag2
5 sdiag1 = sdiag2
sdiag1 = sdiag2
end do
6 diag = diag1-sdiag1*temp(nm1)
yp(n) = (slppn-dx1-sdiag1*yp(nm1))/diag
c
c perform back substitution
c
do 7 i = 2,n
do i = 2,n
ibak = np1-i
7 yp(ibak) = yp(ibak)-temp(ibak)*yp(ibak+1)
yp(ibak) = yp(ibak)-temp(ibak)*yp(ibak+1)
end do
return
c
c too few points
Expand Down Expand Up @@ -681,7 +683,8 @@ subroutine curvp1 (n,x,y,p,yp,temp,sigma,ierr)
diag1 = diag2
sdiag1 = sdiag2
if (n .eq. 2) go to 2
do 1 i = 2,nm1

do i = 2,nm1
npi = n+i
delx2 = x(i+1)-x(i)
if (delx2 .le. 0.) go to 8
Expand All @@ -693,29 +696,35 @@ subroutine curvp1 (n,x,y,p,yp,temp,sigma,ierr)
temp(i) = sdiag2/diag
dx1 = dx2
diag1 = diag2
1 sdiag1 = sdiag2
sdiag1 = sdiag2
end do

2 delx2 = p-(x(n)-x(1))
dx2 = (y(1)-y(n))/delx2
call terms (diag2,sdiag2,sigmap,delx2)
yp(n) = dx2-dx1
temp(nm1) = temp(2*n-1)-temp(nm1)

if (n .eq. 2) go to 4
c
c perform first step of back substitution
c
do 3 i = 3,n
do i = 3,n
ibak = np1-i
npibak =n+ibak
yp(ibak) = yp(ibak)-temp(ibak)*yp(ibak+1)
3 temp(ibak) =temp(npibak)-temp(ibak)*temp(ibak+1)
temp(ibak) =temp(npibak)-temp(ibak)*temp(ibak+1)
end do

4 yp(n) = (yp(n)-sdiag2*yp(1)-sdiag1*yp(nm1))/
* (diag1+diag2+sdiag2*temp(1)+sdiag1*temp(nm1))
c
c perform second step of back substitution
c
ypn = yp(n)
do 5 i = 1,nm1
5 yp(i) = yp(i)+temp(i)*ypn
do i = 1,nm1
yp(i) = yp(i)+temp(i)*ypn
end do
return
c
c too few points
Expand All @@ -733,6 +742,8 @@ subroutine curvp1 (n,x,y,p,yp,temp,sigma,ierr)
8 ierr = 3
return
end

!--------------------------------------------------------------------------------
function curvp2 (t,n,x,y,p,yp,sigma)
c
integer:: n
Expand Down Expand Up @@ -3199,6 +3210,9 @@ subroutine snhcsh (sinhm,coshm,x,isw)
if (isw .eq. 3) sinhm = (expx-1./expx)/(ax+ax)-1.
return
end


!--------------------------------------------------------------------------------
subroutine surf1 (m,n,x,y,z,iz,zx1,zxm,zy1,zyn,zxy11,
* zxym1,zxy1n,zxymn,islpsw,zp,temp,
* sigma,ierr)
Expand Down Expand Up @@ -3522,54 +3536,64 @@ subroutine surf1 (m,n,x,y,z,iz,zx1,zxm,zy1,zyn,zxy11,
del1 = x(2)-x(1)
if (del1 .le. 0.) go to 47
deli = 1./del1
do 38 j = 1,n
do j = 1,n
zp(2,j,2) = deli*(z(2,j)-z(1,j))
38 zp(2,j,3) = deli*(zp(2,j,1)-zp(1,j,1))
zp(2,j,3) = deli*(zp(2,j,1)-zp(1,j,1))
end do
call terms (diag1,sdiag1,sigmax,del1)
diagi = 1./diag1
do 39 j = 1,n
do j = 1,n
zp(1,j,2) = diagi*(zp(2,j,2)-zp(1,j,2))
39 zp(1,j,3) = diagi*(zp(2,j,3)-zp(1,j,3))
zp(1,j,3) = diagi*(zp(2,j,3)-zp(1,j,3))
end do

temp(n+1) = diagi*sdiag1
if (m .eq. 2) go to 43
do 42 i = 2,mm1
do i = 2,mm1
im1 = i-1
ip1 = i+1
npi = n+i
del2 = x(ip1)-x(i)
if (del2 .le. 0.) go to 47
deli = 1./del2
do 40 j = 1,n
do j = 1,n
zp(ip1,j,2) = deli*(z(ip1,j)-z(i,j))
40 zp(ip1,j,3) = deli*(zp(ip1,j,1)-zp(i,j,1))
zp(ip1,j,3) = deli*(zp(ip1,j,1)-zp(i,j,1))
end do
call terms (diag2,sdiag2,sigmax,del2)
diagin = 1./(diag1+diag2-sdiag1*temp(npi-1))
do 41 j = 1,n
do j = 1,n
zp(i,j,2) = diagin*(zp(ip1,j,2)-zp(i,j,2)-
* sdiag1*zp(im1,j,2))
41 zp(i,j,3) = diagin*(zp(ip1,j,3)-zp(i,j,3)-
* sdiag1*zp(im1,j,3))
zp(i,j,3) = diagin*(zp(ip1,j,3)-zp(i,j,3)-
* sdiag1*zp(im1,j,3))
end do
temp(npi) = diagin*sdiag2
diag1 = diag2
42 sdiag1 = sdiag2
sdiag1 = sdiag2
end do

43 diagin = 1./(diag1-sdiag1*temp(npm-1))
do 44 j = 1,n
do j = 1,n
npmpj = npm+j
zp(m,j,2) = diagin*(temp(npmpj)-zp(m,j,2)-
* sdiag1*zp(mm1,j,2))
44 zp(m,j,3) = diagin*(temp(j)-zp(m,j,3)-
* sdiag1*zp(mm1,j,3))
zp(m,j,3) = diagin*(temp(j)-zp(m,j,3)-
* sdiag1*zp(mm1,j,3))
end do
c
c perform back substitution
c
do 45 i = 2,m
do i = 2,m
ibak = mp1-i
ibakp1 = ibak+1
npibak = n+ibak
t = temp(npibak)
do 45 j = 1,n
do j = 1,n
zp(ibak,j,2) = zp(ibak,j,2)-t*zp(ibakp1,j,2)
45 zp(ibak,j,3) = zp(ibak,j,3)-t*zp(ibakp1,j,3)
zp(ibak,j,3) = zp(ibak,j,3)-t*zp(ibakp1,j,3)
end do
end do
return
c
c too few points
Expand All @@ -3582,6 +3606,8 @@ subroutine surf1 (m,n,x,y,z,iz,zx1,zxm,zy1,zyn,zxy11,
47 ierr = 2
return
end

!--------------------------------------------------------------------------------
function surf2 (xx,yy,m,n,x,y,z,iz,zp,sigma)
c
integer:: m,n,iz
Expand Down
23 changes: 16 additions & 7 deletions src/aorsa2dMain.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -866,9 +866,9 @@ program aorsa2dMain
real :: rho_thresh


complex :: aij, beta
! complex :: aij, beta
real, dimension(:,:), allocatable :: work
complex, dimension(:,:), allocatable :: workc
! complex, dimension(:,:), allocatable :: workc

complex, dimension(ndfmax) :: row,rowk

Expand Down Expand Up @@ -911,7 +911,7 @@ program aorsa2dMain
& yy(mkdim1 : mkdim2, 1 : nymx) )

allocate (work(nxmx, nymx))
allocate (workc(nxmx, nymx))
! allocate (workc(nxmx, nymx))

allocate ( psi(nxmx, nymx), rho(nxmx, nymx), theta(nxmx, nymx),
& rhohatx(nxmx, nymx), rhohaty(nxmx, nymx), theta0(nxmx, nymx),
Expand Down Expand Up @@ -2255,6 +2255,13 @@ program aorsa2dMain
uzy(i, j) = byn(i, j)
uzz(i, j) = bzn(i, j)

! Rotation matrix
!
! 1-bxn^2
! U =
! bxn byn bzn
!

dxuxx(i,j) = dxsqx
dxxuxx(i,j) = dxxsqx
dyuxx(i,j) = dysqx
Expand Down Expand Up @@ -3321,7 +3328,8 @@ program aorsa2dMain
xjz(i,j) = xjantz * shapey * gausspsi * Iz

if (myid==0) then
write(*,*) "gausspsi", i,j,gausspsi,xjz(i,j)
write(*,*) "gausspsi", i,j,gausspsi,xjz(i,j),psiant,
& psi(i,j), dpsiant,xjantz,shapey,Iz,deltay,yant_max
end if
end if

Expand Down Expand Up @@ -3376,7 +3384,7 @@ program aorsa2dMain
end if

! Set k_phi = n_phi/R_ant where Rant is the location of the current peak
write(*,*) 'current location',i_cur_max,j_cur_max,capr(i_cur_max),
write(*,*) 'current location',i_cur_max,j_cur_max,
& cur_mod_max,capr(i_cur_max),y(j_cur_max)
xkphia = nphi / capr(i_cur_max)
npar = xkphia / xk0
Expand Down Expand Up @@ -6568,7 +6576,7 @@ program aorsa2dMain
p_amat(ipos + 1) = fgk
p_amat(ipos + 2) = frk
if (icnc1 .eq. 1) then
if (icnc1 .eq. 1) then !if first global row
p_brhs(ipos) = xb(i,j) !xb,c,d = current source
p_brhs(ipos + 1) = xc(i,j)
p_brhs(ipos + 2) = xd(i,j)
Expand Down Expand Up @@ -7263,6 +7271,7 @@ program aorsa2dMain
end do
end do
!JCW Why is this being done, just for plots/ doesnt work for mirror?
*---------------------------------------------------
* Interpolate back onto (R, Z) grid for plotting
*---------------------------------------------------
Expand Down Expand Up @@ -11702,7 +11711,7 @@ program aorsa2dMain
deallocate (xkperp2_slow, xkperp2_fast)
deallocate (xkprl_a, P_a)
deallocate (work)
deallocate (workc)
! deallocate (workc)
!if not static
deallocate (workn)
deallocate (rwork)
Expand Down
12 changes: 6 additions & 6 deletions src/eqdsk_plot.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -784,14 +784,14 @@ subroutine eqdsk_plot
nxmx, nymx, nlevmax, title, titx, tity, scalex)


title = 'bmod_mid surfaces'
title = 'bmod_{mid} surfaces'
call ezconc(capr, y, bmod_mid, ff, nnodex, nnodey, numb, &
nxmx, nymx, nlevmax, title, titx, tity, iflag, scalex)
if (iflag .eq. 0) call boundary(capr, y, rho, ff, nnodex, &
nnodey, numb, nxmx, nymx, nlevmax, title, titx, tity,scalex)


title = 'capr_{bpol-mid2} surfaces'
title = '(R*Bpol)_{mid2} surfaces'
call ezconc(capr, y, capr_bpol_mid2, ff, nnodex, nnodey, numb, &
nxmx, nymx, nlevmax, title, titx, tity, iflag, scalex)
if (iflag .eq. 0) call boundary(capr, y, rho, ff, nnodex, &
Expand Down Expand Up @@ -838,15 +838,15 @@ subroutine eqdsk_plot
nnodex, nxmx)


title= 'Flux average bmod_mid'
title= 'Flux average bmod_{mid}'
titll= 'bmod_mid (T)'
titlr=' '
call ezplot1(title, titll, titlr, rhon_half, bmod_midavg, &
nnoderho_half, nrhomax)


title= 'Flux average capr_bpol'
titll= 'r*bpol (mT)'
title= 'Flux average R*B_{pol}'
titll= '<r*bpol> (m-T)'
titlr=' '
call ezplot1(title, titll, titlr, rhon_half, capr_bpol_midavg, &
nnoderho_half, nrhomax)
Expand Down Expand Up @@ -894,7 +894,7 @@ subroutine eqdsk_plot



title = 'grad_parallel B'
title = 'grad_{||} B'
call ezconc(capr, y, gradprlb, ff, nnodex, nnodey, numb, &
nxmx, nymx, nlevmax, title, titx, tity,iflag, scalex)
if (iflag .eq. 0) call boundary (capr, y, rho, ff, nnodex, &
Expand Down
Loading

0 comments on commit 641606b

Please sign in to comment.