Skip to content

Commit 4568a5c

Browse files
d-burgd-burg
committed
WIP: looks roughly correct but yields wrong answers
1 parent dce3e95 commit 4568a5c

File tree

1 file changed

+62
-45
lines changed

1 file changed

+62
-45
lines changed

‎stride/sym.f‎

Lines changed: 62 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -34,8 +34,8 @@ SUBROUTINE symmetrize(dp)
3434
REAL(r8) :: resnum,resm_sing,min_diff,dtheta
3535
INTEGER, DIMENSION(msing) :: resm
3636
LOGICAL :: is_hermitian
37-
REAL(r8) :: tolerance, dx
38-
REAL(r8), DIMENSION(0:mtheta) :: r,z,theta
37+
REAL(r8) :: tolerance,dx,angle,chi1
38+
REAL(r8), DIMENSION(:), ALLOCATABLE :: r,z,theta
3939
REAL(r8), DIMENSION(sq%mx+1) :: I_psin
4040
REAL(r8),DIMENSION(msing) :: L,f_L,f_S,nu_L,nu_S,DI,J,rho,shear
4141

@@ -51,14 +51,19 @@ SUBROUTINE symmetrize(dp)
5151
REAL(r8), DIMENSION(2,2) :: w
5252
TYPE(spline_type) :: spl
5353

54+
WRITE(*,*)"mthsurf=",mthsurf
55+
56+
mthsurf=40*mthsurf0*MAX(ABS(mlow),ABS(mhigh))
57+
ALLOCATE(r(0:mthsurf),z(0:mthsurf),theta(0:mthsurf))
58+
5459
npsi = sq%mx
5560

56-
psi_a=0.1 !!!!!!!!!!!!!!
61+
!chi1=twopi*psio
62+
psi_a=psio
63+
WRITE(*,*)"psi_a=",psi_a
5764

5865
WRITE(*,*)"sq%mx=",sq%mx
5966
WRITE(*,*)"npsi=",npsi
60-
!WRITE(*,*)"sq%xs=",sq%xs
61-
6267

6368
! construct PEST3 matching data (keep synced with RDCON!)
6469
A_prime=0.0
@@ -94,26 +99,11 @@ SUBROUTINE symmetrize(dp)
9499
resm(ising)=resm_sing
95100
ENDDO
96101

97-
!WRITE(*,*)"respsi=",respsi
98102
WRITE(*,*)"resnum=",resnum
99103
WRITE(*,*)"resm=",resm
100104

101-
! STILL NEED psi_a
102-
!I = 0
103-
!DO ipsi=1,npsi
104-
! !CALL spline_eval(sq,psi(ipsi),0)
105-
! I = I + 2.0*(sq%f(4) / (sq%fs(ipsi,1)/twopi)) !DB this is F!
106-
!ENDDO
107-
105+
! Calculate I(psi_N) integral
108106
I_psin(1) = 0.0 ! integral at x=0 is 0
109-
110-
!WRITE(*,*)"sq%fs(:,4)=",sq%fs(:,4)
111-
!WRITE(*,*)"SIZE(sq%fs(:,4))=",SIZE(sq%fs(:,4))
112-
113-
!WRITE(*,*)"sq%fs(npsi,4)=",sq%fs(npsi,4)
114-
115-
!WRITE(*,*)"sq%fs(5,4)=",sq%fs(5,4)
116-
117107
DO ipsi = 2, npsi
118108
I_psin(ipsi) = 0.0
119109
DO jpsi = 1, ipsi-1
@@ -123,14 +113,19 @@ SUBROUTINE symmetrize(dp)
123113
END DO
124114
END DO
125115

126-
!WRITE(*,*)"I_psin=",I_psin
116+
! Prepare theta quantities for J(psi_N) integral
117+
CALL spline_alloc(spl,mthsurf,4)
118+
theta=(/(itheta,itheta=0,mthsurf)/)/REAL(mthsurf,r8)
119+
spl%xs=theta
120+
psifac=psilim
127121

122+
! Loop across rational surfaces to evaluate remaining quantities
128123
DO ising=1,msing
129124
!resnum=NINT(sing(ising)%q*nn)-mlow+1
130125
respsi=sing(ising)%psifac
131126
WRITE(*,*)"respsi=",respsi
132127

133-
!!! FIND CORRECT I_PSIN INDEX
128+
! Find correct I(psiN) index
134129
min_diff = abs(sq%xs(1) - respsi)
135130
idx = 1
136131

@@ -143,51 +138,65 @@ SUBROUTINE symmetrize(dp)
143138

144139
WRITE(*,*)"idx=",idx
145140

146-
!CALL spline_eval(sq,respsi,1)
141+
! Evaluate splines on rational surface
142+
CALL spline_eval(sq,respsi,1)
143+
CALL spline_eval(locstab,respsi,1)
144+
145+
WRITE(*,*)"SIZE(theta)=",SIZE(theta)
146+
WRITE(*,*)"mthsurf=",mthsurf
147+
148+
! Calculate J(psi_N) contour integral
147149
dtheta = twopi / real(mthsurf)
150+
DO itheta=0,mthsurf
151+
CALL bicube_eval(rzphi,psilim,theta(itheta),1)
152+
rfac=SQRT(rzphi%f(1))
153+
angle=twopi*(theta(itheta)+rzphi%f(2))
154+
r(itheta)=ro+rfac*COS(angle)
155+
z(itheta)=zo+rfac*SIN(angle)
156+
jac=rzphi%f(4)
148157

149-
!DO itheta=0,mthsurf
150-
! CALL bicube_eval(rzphi,respsi,theta(itheta),1)
151-
! jac=rzphi%f(4)
152-
! w(1,1)=(1+rzphi%fy(2))*twopi**2*rfac*r(itheta)/jac
153-
! w(1,2)=-rzphi%fy(1)*pi*r(itheta)/(rfac*jac)
154-
! delpsi=SQRT(w(1,1)**2+w(1,2)**2)
155-
! WRITE(*,*)"delpsi=",delpsi
158+
w(1,1)=(1+rzphi%fy(2))*twopi**2*rfac*r(itheta)/jac
159+
w(1,2)=-rzphi%fy(1)*pi*r(itheta)/(rfac*jac)
156160

157-
! J(ising)=J(ising)+(1.0/delpsi**2)*0.5*dtheta!jac!*(1.0/twopi) !DB !!!
158-
!ENDDO
161+
delpsi=SQRT(w(1,1)**2+w(1,2)**2)
162+
163+
J(ising)=J(ising)+
164+
$ (1.0/(delpsi**2))*(dtheta/twopi)!/jac !DB !!
165+
ENDDO
159166

160167
WRITE(*,*)"J=",J
161168

162-
rho(ising) = (J(ising) * (sq%fs(idx,1)/twopi))/sq%f(4) !DB
163-
shear(ising)=sing(ising)%q1
169+
rho(ising) = (J(ising) * (sq%f(1)/twopi))/sq%f(4) !DB
170+
shear(ising) = sing(ising)%q1
164171
WRITE(*,*)"rho=",rho
172+
WRITE(*,*)"shear(ising)=",shear(ising)
173+
165174
WRITE(*,*)"sq%f(4)=",sq%f(4)
166175

176+
! Combine integrated quantities into L(psi_N)
167177
L(ising)=I_psin(idx)*(J(ising)*((resm(ising)*
168-
$ (sq%fs(idx,1)/twopi))/sq%f(4))**2+(nn * psi_a)**2 ) !DB
178+
$ (sq%f(1)/twopi))/sq%f(4))**2+(nn * psi_a)**2 ) !DB
169179

170180
WRITE(*,*)"L=",L
171-
WRITE(*,*)"sq%fs(respsi,1)=",sq%fs(respsi,1)
172-
WRITE(*,*)"sq%fs(idx,1)=",sq%fs(idx,1)
173-
174181
! CALL mercier_mod
175-
DI(ising) = locstab%fs(idx,1)/sq%xs(idx) !DB
182+
DI(ising) = locstab%f(1)/respsi !DB
176183
WRITE(*,*)"DI=",DI
177184

178185
nu_L(ising) = 0.5 - SQRT(-DI(ising)) !DB
179186
nu_S(ising) = 0.5 + SQRT(-DI(ising)) !DB
180187

181-
f_L(ising) = (rho(ising) ** nu_L(ising)) * (((nu_S(ising) -
188+
f_L(ising) = (rho(ising) ** nu_L(ising)) * (((nu_S(ising)-
182189
$ nu_L(ising)) /L(ising))**0.5)*shear(ising)*resm(ising) !DB
183-
f_S(ising) = (rho(ising) ** nu_S(ising)) * (((nu_S(ising) -
190+
f_S(ising) = (rho(ising) ** nu_S(ising)) * (((nu_S(ising)-
184191
$ nu_L(ising)) /L(ising))**0.5)*shear(ising)*resm(ising) !DB
185192

193+
! First component of vector*matrix*vector multiply
186194
A_prime_tmp = MATMUL(A_prime,f_L) !DB
187195
B_prime_tmp = MATMUL(B_prime,f_L) !DB
188196
Gamma_prime_tmp = MATMUL(Gamma_prime,f_L) !DB
189197
Delta_prime_tmp = MATMUL(Delta_prime,f_L) !DB
190198

199+
! Second component of vector*matrix*vector multiply
191200
A_prime_sym = MATMUL(RESHAPE([f_S], [1,msing]),
192201
$ RESHAPE([A_prime_tmp], [1,msing])) !DB
193202
B_prime_sym = MATMUL(RESHAPE([f_S], [1,msing]),
@@ -198,7 +207,15 @@ SUBROUTINE symmetrize(dp)
198207
$ RESHAPE([Delta_prime_tmp], [1,msing])) !DB
199208
ENDDO
200209

201-
tolerance = 1.0e-2
210+
WRITE(*,*)"Delta_prime=",Delta_prime
211+
WRITE(*,*)"Delta_prime_sym=",Delta_prime_sym
212+
WRITE(*,*)"Delta_diff=",abs(Delta_prime_sym -
213+
$ transpose(conjg(Delta_prime_sym)))
214+
!WRITE(*,*)"B_prime=",B_prime
215+
!WRITE(*,*)"Gamma_prime_sym=",Gamma_prime_sym
216+
WRITE(*,*)"B_Gamma_diff=",(B_prime_sym -
217+
$ transpose(conjg(Gamma_prime_sym)))
218+
tolerance = 1.0e-01
202219

203220
is_hermitian = all(abs(A_prime_sym -
204221
$ transpose(conjg(A_prime_sym))) < tolerance)
@@ -207,14 +224,14 @@ SUBROUTINE symmetrize(dp)
207224
ELSE
208225
WRITE(*,*), "A_prime is not Hermitian"
209226
END IF
210-
is_hermitian = all(abs(B_prime_sym -
227+
is_hermitian = all(abs(Gamma_prime_sym -
211228
$ transpose(conjg(B_prime_sym))) < tolerance)
212229
IF (is_hermitian) then
213230
WRITE(*,*), "B_prime is Hermitian"
214231
ELSE
215232
WRITE(*,*), "B_prime is not Hermitian"
216233
END IF
217-
is_hermitian = all(abs(Gamma_prime_sym -
234+
is_hermitian = all(abs(B_prime_sym -
218235
$ transpose(conjg(Gamma_prime_sym))) < tolerance)
219236
IF (is_hermitian) then
220237
WRITE(*,*), "Gamma_prime_sym is Hermitian"

0 commit comments

Comments
 (0)