Skip to content

Commit

Permalink
add dbgout for MHD forces
Browse files Browse the repository at this point in the history
  • Loading branch information
@jons-pf
jons-pf committed Oct 15, 2024
1 parent 70e311b commit d550d20
Showing 1 changed file with 148 additions and 22 deletions.
170 changes: 148 additions & 22 deletions Sources/General/forces.f
View file
Original file line number Diff line number Diff line change
@@ -1,13 +1,15 @@
SUBROUTINE forces_par
USE vmec_main, p5 => cp5
USE realspace
USE vforces, ru12 => pazmn_e, zu12 => parmn_e,
USE vforces, ru12 => pazmn_e, zu12 => parmn_e,
& pazmn_e => pazmn_e, parmn_e => parmn_e,
& lv_e => pcrmn_e, lu_e => pczmn_e, lu_o => pczmn_o,
& pcrmn_e => pcrmn_e, pczmn_e => pczmn_e,
& pczmn_o => pczmn_o
USE parallel_include_module
USE timer_sub
USE dbgout

IMPLICIT NONE
C-----------------------------------------------
C L o c a l P a r a m e t e r s
Expand All @@ -18,14 +20,71 @@ SUBROUTINE forces_par
C-----------------------------------------------
INTEGER :: l, ndim
INTEGER :: i, j, k, nsmin, nsmax
REAL(dp), DIMENSION(:,:), POINTER ::
REAL(dp), DIMENSION(:,:), POINTER ::
& bsqr, gvvs, guvs, guus
REAL(dp), ALLOCATABLE, DIMENSION(:) :: bcastbuf
LOGICAL :: dbg_forces
C-----------------------------------------------
IF (.NOT.lactive .AND. .NOT.lfreeb) RETURN

CALL second0 (tforon)

dbg_forces = open_dbg_context("forces", num_eqsolve_retries)
if (dbg_forces) then

call add_real_3d("lu_e", ns, nzeta, ntheta3, lu_e, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("lv_e_in", ns, nzeta, ntheta3, lv_e, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("guu_in", ns, nzeta, ntheta3, pguu, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("guv_in", ns, nzeta, ntheta3, pguv, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("gvv_in", ns, nzeta, ntheta3, pgvv, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("ru12", ns, nzeta, ntheta3, ru12, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("zu12", ns, nzeta, ntheta3, zu12, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("brmn_e_in", ns, nzeta, ntheta3, pbrmn_e, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("bzmn_e_in", ns, nzeta, ntheta3, pbzmn_e, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("ru0", ns, nzeta, ntheta3, pru0, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("zu0", ns, nzeta, ntheta3, pzu0, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("rcon0", ns, nzeta, ntheta3, prcon0, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("zcon0", ns, nzeta, ntheta3, pzcon0, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("gcon", ns, nzeta, ntheta3, pgcon, &
& order = (/ 2, 3, 1 /) )

call add_real_4d("r1", ns, 2, nzeta, ntheta3, pr1, &
& order=(/ 3, 4, 1, 2 /) )
call add_real_4d("z1", ns, 2, nzeta, ntheta3, pz1, &
& order=(/ 3, 4, 1, 2 /) )
call add_real_4d("ru", ns, 2, nzeta, ntheta3, pru, &
& order=(/ 3, 4, 1, 2 /) )
call add_real_4d("zu", ns, 2, nzeta, ntheta3, pzu, &
& order=(/ 3, 4, 1, 2 /) )
call add_real_4d("rcon_in", ns, 2, nzeta, ntheta3, prcon, &
& order=(/ 3, 4, 1, 2 /) )
call add_real_4d("zcon_in", ns, 2, nzeta, ntheta3, pzcon, &
& order=(/ 3, 4, 1, 2 /) )

if (lthreed) then
call add_real_4d("rv", ns, 2, nzeta, ntheta3, prv, &
& order=(/ 3, 4, 1, 2 /) )
call add_real_4d("zv", ns, 2, nzeta, ntheta3, pzv, &
& order=(/ 3, 4, 1, 2 /) )
else
call add_null("rv")
call add_null("zv")
end if
end if ! dump_forces

ndim = 1+nrzt

nsmin=tlglob; nsmax=t1rglob
Expand All @@ -37,7 +96,7 @@ SUBROUTINE forces_par
lu_e(:,1) = 0; lv_e(:,1) = 0
pguu(:,1) = 0; pguv(:,1) = 0; pgvv(:,1) = 0

DO l = nsmin, nsmax
DO l = nsmin, nsmax
guus(:,l) = pguu(:,l)*pshalf(:,l)
guvs(:,l) = pguv(:,l)*pshalf(:,l)
gvvs(:,l) = pgvv(:,l)* pshalf(:,l)
Expand Down Expand Up @@ -86,10 +145,10 @@ SUBROUTINE forces_par
pbzmn_e(:,l) = p5*(pbzmn_e(:,l) + pbzmn_e(:,l+1))
END DO

parmn_e(:,ns) = - parmn_e(:,ns) + p5*lv_e(:,ns)
parmn_e(:,ns) = - parmn_e(:,ns) + p5*lv_e(:,ns)
pazmn_e(:,ns) = - pazmn_e(:,ns)
pbrmn_e(:,ns) = p5*pbrmn_e(:,ns)
pbzmn_e(:,ns) = p5*pbzmn_e(:,ns)
pbzmn_e(:,ns) = p5*pbzmn_e(:,ns)

nsmin=tlglob; nsmax=t1rglob
DO l = nsmin, nsmax
Expand Down Expand Up @@ -120,9 +179,19 @@ SUBROUTINE forces_par
& + p5*lv_e(:,ns)
pazmn_o(:,ns) = - pazmn_o(:,ns) + pru(:,ns,0)*bsqr(:,ns)
pbrmn_o(:,ns) = p5*pbrmn_o(:,ns)
pbzmn_o(:,ns) = p5*pbzmn_o(:,ns)
pbzmn_o(:,ns) = p5*pbzmn_o(:,ns)
lu_o(:,ns) = lu_o(:,ns)

if (dbg_forces) then
! save data before it gets overwritten again
! due to array re-usage

call add_real_3d("lu_o", ns, nzeta, ntheta3, lu_o, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("lv_e_out", ns, nzeta, ntheta3, lv_e, &
& order = (/ 2, 3, 1 /) )
end if

nsmin=tlglob; nsmax=trglob
DO l = nsmin, nsmax
pguu(:,l) = pguu(:,l) * psqrts(:,l)**2
Expand Down Expand Up @@ -207,6 +276,63 @@ SUBROUTINE forces_par
pzcon(:,l,1) = pzcon(:,l,0) * psqrts(:,l)
END DO
#endif

if (dbg_forces) then
call add_real_3d("armn_e", ns, nzeta, ntheta3, parmn_e, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("armn_o", ns, nzeta, ntheta3, parmn_o, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("brmn_e", ns, nzeta, ntheta3, pbrmn_e, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("brmn_o", ns, nzeta, ntheta3, pbrmn_o, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("azmn_e", ns, nzeta, ntheta3, pazmn_e, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("azmn_o", ns, nzeta, ntheta3, pazmn_o, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("bzmn_e", ns, nzeta, ntheta3, pbzmn_e, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("bzmn_o", ns, nzeta, ntheta3, pbzmn_o, &
& order = (/ 2, 3, 1 /) )
if (lthreed) then
call add_real_3d("crmn_e", ns, nzeta, ntheta3, pcrmn_e, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("crmn_o", ns, nzeta, ntheta3, pcrmn_o, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("czmn_e", ns, nzeta, ntheta3, pczmn_e, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("czmn_o", ns, nzeta, ntheta3, pczmn_o, &
& order = (/ 2, 3, 1 /) )
else
call add_null("crmn_e")
call add_null("crmn_o")
call add_null("czmn_e")
call add_null("czmn_o")
end if

call add_real_3d("guu_out", ns, nzeta, ntheta3, pguu, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("guus", ns, nzeta, ntheta3, guus, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("guv_out", ns, nzeta, ntheta3, pguv, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("guvs", ns, nzeta, ntheta3, guvs, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("gvv_out", ns, nzeta, ntheta3, pgvv, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("gvvs", ns, nzeta, ntheta3, gvvs, &
& order = (/ 2, 3, 1 /) )
call add_real_3d("bsqr", ns, nzeta, ntheta3, bsqr, &
& order = (/ 2, 3, 1 /) )

call add_real_4d("rcon_out", ns, 2, nzeta, ntheta3, prcon, &
& order=(/ 3, 4, 1, 2 /) )
call add_real_4d("zcon_out", ns, 2, nzeta, ntheta3, pzcon, &
& order=(/ 3, 4, 1, 2 /) )

call close_dbg_out()
end if

CALL second0 (tforoff)
timer(tfor) = timer(tfor) + (tforoff - tforon)
forces_time = timer(tfor)
Expand All @@ -216,7 +342,7 @@ END SUBROUTINE forces_par
SUBROUTINE forces
USE vmec_main, p5 => cp5
USE realspace
USE vforces, ru12 => azmn_e, zu12 => armn_e,
USE vforces, ru12 => azmn_e, zu12 => armn_e,
& azmn_e => azmn_e, armn_e => armn_e,
& lv_e => crmn_e, lu_e => czmn_e, lu_o => czmn_o,
& crmn_e => crmn_e, czmn_e => czmn_e, czmn_o => czmn_o
Expand All @@ -231,7 +357,7 @@ SUBROUTINE forces
C L o c a l V a r i a b l e s
C-----------------------------------------------
INTEGER :: l, ndim
REAL(dp), DIMENSION(:), POINTER ::
REAL(dp), DIMENSION(:), POINTER ::
& bsqr, gvvs, guvs, guus

C-----------------------------------------------
Expand Down Expand Up @@ -281,7 +407,7 @@ SUBROUTINE forces
azmn_o(1:ndim) = azmn_e(1:ndim) *shalf
brmn_o(1:ndim) = brmn_e(1:ndim) *shalf
bzmn_o(1:ndim) = bzmn_e(1:ndim) *shalf

!
! CONSTRUCT CYLINDRICAL FORCE KERNELS
! NOTE: presg(ns+1) == 0, AND WILL BE "FILLED IN" AT EDGE
Expand Down Expand Up @@ -341,25 +467,25 @@ SUBROUTINE forces
guvs(l) = p5*(guvs(l) + guvs(l+1))
END DO

brmn_e(:nrzt) = brmn_e(:nrzt)
brmn_e(:nrzt) = brmn_e(:nrzt)
& - (guv(:nrzt)*rv(:nrzt,0) + guvs(:nrzt)*rv(:nrzt,1))
bzmn_e(:nrzt) = bzmn_e(:nrzt)
bzmn_e(:nrzt) = bzmn_e(:nrzt)
& - (guv(:nrzt)*zv(:nrzt,0) + guvs(:nrzt)*zv(:nrzt,1))
crmn_e(:nrzt) = guv(:nrzt) *ru(:nrzt,0)
crmn_e(:nrzt) = guv(:nrzt) *ru(:nrzt,0)
& + gvv(:nrzt) *rv(:nrzt,0)
& + gvvs(:nrzt)*rv(:nrzt,1) + guvs(:nrzt)*ru(:nrzt,1)
czmn_e(:nrzt) = guv(:nrzt) *zu(:nrzt,0)
czmn_e(:nrzt) = guv(:nrzt) *zu(:nrzt,0)
& + gvv(:nrzt) *zv(:nrzt,0)
& + gvvs(:nrzt)*zv(:nrzt,1) + guvs(:nrzt)*zu(:nrzt,1)
guv(:nrzt) = guv(:nrzt) *sqrts(:nrzt)*sqrts(:nrzt)
brmn_o(:nrzt) = brmn_o(:nrzt)
brmn_o(:nrzt) = brmn_o(:nrzt)
& - (guvs(:nrzt)*rv(:nrzt,0) + guv(:nrzt)*rv(:nrzt,1))
bzmn_o(:nrzt) = bzmn_o(:nrzt)
bzmn_o(:nrzt) = bzmn_o(:nrzt)
& - (guvs(:nrzt)*zv(:nrzt,0) + guv(:nrzt)*zv(:nrzt,1))
crmn_o(:nrzt) = guvs(:nrzt)*ru(:nrzt,0)
crmn_o(:nrzt) = guvs(:nrzt)*ru(:nrzt,0)
& + gvvs(:nrzt)*rv(:nrzt,0)
& + bsqr(:nrzt)*rv(:nrzt,1) + guv(:nrzt) *ru(:nrzt,1)
czmn_o(:nrzt) = guvs(:nrzt)*zu(:nrzt,0)
czmn_o(:nrzt) = guvs(:nrzt)*zu(:nrzt,0)
& + gvvs(:nrzt)*zv(:nrzt,0)
& + bsqr(:nrzt)*zv(:nrzt,1) + guv(:nrzt) *zu(:nrzt,1)
ENDIF
Expand All @@ -368,13 +494,13 @@ SUBROUTINE forces
! ASSIGN EDGE FORCES (JS = NS) FOR FREE BOUNDARY CALCULATION
!
IF (ivac .ge. 1) THEN
armn_e(ns:nrzt:ns) = armn_e(ns:nrzt:ns)
armn_e(ns:nrzt:ns) = armn_e(ns:nrzt:ns)
& + zu0(ns:nrzt:ns)*rbsq(1:nznt)
armn_o(ns:nrzt:ns) = armn_o(ns:nrzt:ns)
armn_o(ns:nrzt:ns) = armn_o(ns:nrzt:ns)
& + zu0(ns:nrzt:ns)*rbsq(1:nznt)
azmn_e(ns:nrzt:ns) = azmn_e(ns:nrzt:ns)
azmn_e(ns:nrzt:ns) = azmn_e(ns:nrzt:ns)
& - ru0(ns:nrzt:ns)*rbsq(1:nznt)
azmn_o(ns:nrzt:ns) = azmn_o(ns:nrzt:ns)
azmn_o(ns:nrzt:ns) = azmn_o(ns:nrzt:ns)
& - ru0(ns:nrzt:ns)*rbsq(1:nznt)
! fz00_edge = SUM(wint(ns:nrzt:ns)*ru0(ns:nrzt:ns)*rbsq(1:nznt))
ENDIF
Expand All @@ -397,7 +523,7 @@ SUBROUTINE forces
rcon(:nrzt,1) = rcon(:nrzt,0) * sqrts(:nrzt)
zcon(:nrzt,1) = zcon(:nrzt,0) * sqrts(:nrzt)
#endif

CALL second0 (tforoff)
timer(tfor) = timer(tfor) + (tforoff - tforon)

Expand Down

0 comments on commit d550d20

Please sign in to comment.