Skip to content

Commit

Permalink
some random debugging snapshot?
Browse files Browse the repository at this point in the history
  • Loading branch information
jons-pf committed Jul 29, 2024
1 parent 49a47d3 commit 37feb94
Show file tree
Hide file tree
Showing 4 changed files with 43 additions and 28 deletions.
58 changes: 34 additions & 24 deletions src/NESTOR/analyt.f90
Original file line number Diff line number Diff line change
Expand Up @@ -11,13 +11,13 @@
!> @param n_map
!> @param grpmn_m_map
!> @param grpmn_n_map
SUBROUTINE analyt(grpmn, bvec, ivacskip, lasym, m_map, n_map, grpmn_m_map, grpmn_n_map)
SUBROUTINE analyt(grpmn, bvec, ivacskip, lasym, m_map, n_map, grpmn_m_map, grpmn_n_map, ivac)
USE vacmod, vm_grpmn => grpmn
use dbgout
use vmec_main, only: num_eqsolve_retries
IMPLICIT NONE

INTEGER, INTENT(in) :: ivacskip
INTEGER, INTENT(in) :: ivacskip, ivac
REAL(rprec), INTENT(out) :: grpmn(nuv2*mnpd2)
REAL(rprec), INTENT(out) :: bvec(mnpd2)
real(rprec), intent(out) :: m_map(mnpd2)
Expand Down Expand Up @@ -70,16 +70,16 @@ SUBROUTINE analyt(grpmn, bvec, ivacskip, lasym, m_map, n_map, grpmn_m_map, grpmn
! R0P(M): Coefficient of l*T(l-1)+(-) in eq (A17)
! RA1P(M):Coefficient of Tl+(-) in eq (A17)

adp = guu_b + guv_b + gvv_b
adm = guu_b - guv_b + gvv_b
cma = gvv_b - guu_b
sqrtc = two*SQRT(gvv_b)
sqrta = two*SQRT(guu_b)
sqad1u = SQRT(adp)
sqad2u = SQRT(adm)
adp = guu_b + guv_b + gvv_b ! ok
adm = guu_b - guv_b + gvv_b ! ok
cma = gvv_b - guu_b ! ok
sqrtc = two*SQRT(gvv_b) ! ok
sqrta = two*SQRT(guu_b) ! ok
sqad1u = SQRT(adp) ! ok
sqad2u = SQRT(adm) ! ok

! INITIALIZE VECTORS
bvec = 0
bvec = 0 ! ok
IF (ivacskip .eq. 0) THEN
grpmn = 0

Expand All @@ -100,20 +100,20 @@ SUBROUTINE analyt(grpmn, bvec, ivacskip, lasym, m_map, n_map, grpmn_m_map, grpmn
! TLP(M): TL+(-)
! TLP(M)1:T(L-1)+(-)
! TLP(M)2:T(L-2)+(-)
tlp1 = 0
tlm1 = 0
tlp = one/sqad1u*log((sqad1u*sqrtc + adp + cma)/(sqad1u*sqrta - adp + cma))
tlm = one/sqad2u*log((sqad2u*sqrtc + adm + cma)/(sqad2u*sqrta - adm + cma))
tlpm = tlp + tlm
tlp1 = 0 ! ok
tlm1 = 0 ! ok
tlp = one/sqad1u*log((sqad1u*sqrtc + adp + cma)/(sqad1u*sqrta - adp + cma)) ! ok
tlm = one/sqad2u*log((sqad2u*sqrtc + adm + cma)/(sqad2u*sqrta - adm + cma)) ! ok
tlpm = tlp + tlm ! ok

! BEGIN L-SUM IN EQ (A14) TO COMPUTE Imn (and Kmn) INTEGRALS
! NOTE THAT IN THE LOOP OVER L BELOW: L == |m - n| + 2L_A14
! THUS, L BELOW IS THE INDEX OF THE T+- (S+-)
sign1 = 1
fl1 = 0
sign1 = 1 ! ok
fl1 = 0 ! ok

LLOOP: DO l = 0, mf + nf
fl = fl1
fl = fl1 ! ok

! here, tlp/m are available for the current value of l
! --> save into matrix for debugging
Expand All @@ -133,6 +133,12 @@ SUBROUTINE analyt(grpmn, bvec, ivacskip, lasym, m_map, n_map, grpmn_m_map, grpmn
all_slm(l,:) = slm
ENDIF

! if (ivac .eq. 0) then
! ! print *, tlpm
! ! print *, tlp
! ! print *, tlm
! end if

! BEGIN MODE NUMBER (m,n) LOOP
DO n = 0, nf
DO m = 0, mf
Expand All @@ -151,17 +157,21 @@ SUBROUTINE analyt(grpmn, bvec, ivacskip, lasym, m_map, n_map, grpmn_m_map, grpmn
END DO
END DO

! if (ivac .eq. 0) then
! print *, bvec
! end if

! UPDATE "TL's" (FOR L -> L+1) USING EQ (A15)
tlp2 = tlp1
tlm2 = tlm1
tlp1 = tlp
tlm1 = tlm
tlp2 = tlp1 ! ok
tlm2 = tlm1 ! ok
tlp1 = tlp ! ok
tlm1 = tlm ! ok

! next l
fl1 = fl1 + 1
fl1 = fl1 + 1 ! ok

! (-1)**l (next l now)
sign1 = -sign1
sign1 = -sign1 ! ok

tlp = ((sqrtc + sign1*sqrta) - (two*fl1 - one)*cma*tlp1 - fl*adm*tlp2)/(adp*fl1)
tlm = ((sqrtc + sign1*sqrta) - (two*fl1 - one)*cma*tlm1 - fl*adp*tlm2)/(adm*fl1)
Expand Down
5 changes: 5 additions & 0 deletions src/NESTOR/precal.f90
Original file line number Diff line number Diff line change
Expand Up @@ -236,6 +236,11 @@ SUBROUTINE precal

if (open_dbg_context("vac1n_precal", num_eqsolve_retries)) then

call add_real_2d("cosu1", nuv2, mf+1, cosu1)
call add_real_2d("sinu1", nuv2, mf+1, sinu1)
call add_real_2d("cosv1", nuv2, nf+1, cosv1)
call add_real_2d("sinv1", nuv2, nf+1, sinv1)

call add_int("nvper", nvper)
call add_int("nuv_tan", nuv_tan)

Expand Down
6 changes: 3 additions & 3 deletions src/NESTOR/scalpot.f90
Original file line number Diff line number Diff line change
Expand Up @@ -10,15 +10,15 @@
!> @param lasym
!> @param m_map
!> @param n_map
SUBROUTINE scalpot(bvec, amatrix, wint, ivacskip, lasym, m_map, n_map)
SUBROUTINE scalpot(bvec, amatrix, wint, ivacskip, lasym, m_map, n_map, ivac)
USE vacmod, vm_amatrix => amatrix

use dbgout
use vmec_main, only: num_eqsolve_retries

IMPLICIT NONE

INTEGER, INTENT(in) :: ivacskip
INTEGER, INTENT(in) :: ivacskip, ivac
REAL(rprec), INTENT(out) :: bvec(mnpd2), amatrix(mnpd2*mnpd2), m_map(mnpd2), n_map(mnpd2)
REAL(rprec), dimension(nuv2), INTENT(in) :: wint
logical, intent(in) :: lasym
Expand All @@ -34,7 +34,7 @@ SUBROUTINE scalpot(bvec, amatrix, wint, ivacskip, lasym, m_map, n_map)
! BVEC CONTAINS THE TRANSFORM OF THE ANALYTIC SOURCE AND
! GRPMN CONTAINS THE TRANSFORM OF THE NORMAL DERIVATIVE OF THE GREENS FUNCTION [PKM, EQ.(2.15)]
! GReen's function Primed (normal derivative...) and Fourier-transformed to MN mode numbers --> "GR P MN"
CALL analyt (grpmn, bvec, ivacskip, lasym, m_map, n_map, grpmn_m_map_wrt, grpmn_n_map_wrt)
CALL analyt (grpmn, bvec, ivacskip, lasym, m_map, n_map, grpmn_m_map_wrt, grpmn_n_map_wrt, ivac)

IF (ivacskip .ne. 0) THEN
! FOR ivacskip != 0, USE PREVIOUSLY COMPUTED bvecsav FOR SPEED
Expand Down
2 changes: 1 addition & 1 deletion src/NESTOR/vacuum.f90
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,7 @@ SUBROUTINE vacuum(rmnc, rmns, zmns, zmnc, xm, xn, &
! NOTE: all fine up to here against NESTOR.py

! Determine scalar magnetic potential POTVAC
CALL scalpot (potvac, amatrix, wint, ivac_skip, lasym, m_map_wrt, n_map_wrt)
CALL scalpot (potvac, amatrix, wint, ivac_skip, lasym, m_map_wrt, n_map_wrt, ivac)

! stand-alone for debugging: working on scalpot at the moment
! return
Expand Down

0 comments on commit 37feb94

Please sign in to comment.