Skip to content

Carry B_s to the wout file on the full radial mesh #25

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
CharlesCNorton wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Carry B_s to the wout file on the full radial mesh #25

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c24881a
Carry B_s to the wout file on the full radial mesh
CharlesCNorton Sep 3, 2026
5d2cc78
Extrapolate B_s to the boundary from the outermost half point
CharlesCNorton Sep 3, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 10 additions & 4 deletions src/jxbforce.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,8 @@
!> @param bsupv contravariant component of magnetic field \f$B^\zeta\f$
!> @param bsubu covariant component of magnetic field \f$B_\theta\f$
!> @param bsubv covariant component of magnetic field \f$B_\zeta\f$
!> @param bsubsh covariant component of magnetic field \f$B_s\f$ on half grid from bss()
!> @param bsubsh covariant component of magnetic field \f$B_s\f$: on the half grid
!> from bss() on entry, on the full grid on exit
!> @param bsubsu tangential derivate of covariant component of magnetic field \f$\partial B_s / \partial \theta\f$ (?)
!> @param bsubsv tangential derivate of covariant component of magnetic field \f$\partial B_s / \partial \zeta\f$ (?)
!> @param gsqrt Jacobian \f$\sqrt{g}\f$
Expand All @@ -33,7 +34,7 @@ SUBROUTINE jxbforce(bsupu, bsupv, bsubu, bsubv, bsubsh, &
REAL(rprec), DIMENSION(ns,nznt), INTENT(in) :: bsupv
REAL(rprec), DIMENSION(ns,nznt,0:1), TARGET, INTENT(inout) :: bsubu
REAL(rprec), DIMENSION(ns,nznt,0:1), TARGET, INTENT(inout) :: bsubv
REAL(rprec), DIMENSION(ns,nznt), INTENT(in) :: bsubsh
REAL(rprec), DIMENSION(ns,nznt), INTENT(inout) :: bsubsh
REAL(rprec), DIMENSION(ns,nznt,0:1) :: bsubsu
REAL(rprec), DIMENSION(ns,nznt,0:1) :: bsubsv
REAL(rprec), DIMENSION(ns,nznt), INTENT(in) :: gsqrt
Expand Down Expand Up @@ -518,8 +519,13 @@ SUBROUTINE jxbforce(bsupu, bsupv, bsubu, bsubv, bsubsh, &

! Compute end point values for bsubs
bsubs(1,:) = 2*bsubs(2,:) - bsubs(3,:)
!bsubs(ns,:) = 2*bsubs(ns,:) - bsubs(ns-1,:) ! TODO: from ns, ns-1 to ns ???
bsubs(ns,:) = 2*bsubs(ns-1,:) - bsubs(ns-2,:)
! The outermost half point sits between bsubs(ns-1) and the boundary, so
! the boundary value follows from it and one full-mesh point.
bsubs(ns,:) = 2*bsubsh(ns,:) - bsubs(ns-1,:)

! Return the full-mesh B_s to the caller: wrout transforms this array into
! bsubsmns, which the wout file declares on the full mesh.
bsubsh = bsubs



Expand Down
14 changes: 6 additions & 8 deletions src/read_wout_mod.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -727,8 +727,8 @@ SUBROUTINE Compute_Currents(bsubsmnc_i, bsubsmns_i,
! Computes current harmonics for currXmn == sqrt(g)*JsupX, X = u,v
! [Corrected above "JsubX" to "JsupX", JDH 2010-08-16]

! NOTE: bsub(s,u,v)mn are on HALF radial grid
! (in earlier versions, bsubsmn was on FULL radial grid)
! NOTE: bsub(u,v)mn are on the HALF radial grid, bsubsmn on the FULL one,
! as in ORNL-Fusion/LIBSTELL read_wout_mod.f90

!
ohs = (ns_i-1)
Expand All @@ -744,16 +744,15 @@ SUBROUTINE Compute_Currents(bsubsmnc_i, bsubsmns_i,

DO js = 2, ns_i-1
WHERE (MOD(INT(xm_nyq_i),2) .EQ. 1)
t1 = 0.5_dp*(shalf(js+1)*bsubsmns_i(:,js+1) +
& shalf(js) *bsubsmns_i(:,js)) /sfull(js)
t1 = bsubsmns_i(:,js)
bu0 = bsubumnc_i(:,js )/shalf(js)
bu1 = bsubumnc_i(:,js+1)/shalf(js+1)
t2 = ohs*(bu1-bu0)*sfull(js)+0.25_dp*(bu0+bu1)/sfull(js)
bv0 = bsubvmnc_i(:,js )/shalf(js)
bv1 = bsubvmnc_i(:,js+1)/shalf(js+1)
t3 = ohs*(bv1-bv0)*sfull(js)+0.25_dp*(bv0+bv1)/sfull(js)
ELSEWHERE
t1 = 0.5_dp*(bsubsmns_i(:,js+1)+bsubsmns_i(:,js))
t1 = bsubsmns_i(:,js)
t2 = ohs*(bsubumnc_i(:,js+1)-bsubumnc_i(:,js))
t3 = ohs*(bsubvmnc_i(:,js+1)-bsubvmnc_i(:,js))
ENDWHERE
Expand Down Expand Up @@ -783,16 +782,15 @@ SUBROUTINE Compute_Currents(bsubsmnc_i, bsubsmns_i,

DO js = 2, ns_i-1
WHERE (MOD(INT(xm_nyq_i),2) .EQ. 1)
t1 = 0.5_dp*(shalf(js+1)*bsubsmnc_i(:,js+1)
& + shalf(js) *bsubsmnc_i(:,js)) / sfull(js)
t1 = bsubsmnc_i(:,js)
bu0 = bsubumns_i(:,js )/shalf(js+1)
bu1 = bsubumns_i(:,js+1)/shalf(js+1)
t2 = ohs*(bu1-bu0)*sfull(js) + 0.25_dp*(bu0+bu1)/sfull(js)
bv0 = bsubvmns_i(:,js )/shalf(js)
bv1 = bsubvmns_i(:,js+1)/shalf(js+1)
t3 = ohs*(bv1-bv0)*sfull(js)+0.25_dp*(bv0+bv1)/sfull(js)
ELSEWHERE
t1 = 0.5_dp*(bsubsmnc_i(:,js+1) + bsubsmnc_i(:,js))
t1 = bsubsmnc_i(:,js)
t2 = ohs*(bsubumns_i(:,js+1)-bsubumns_i(:,js))
t3 = ohs*(bsubvmns_i(:,js+1)-bsubvmns_i(:,js))
END WHERE
Expand Down
5 changes: 1 addition & 4 deletions src/wrout.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -602,10 +602,7 @@ SUBROUTINE wrout(bsq, gsqrt, bsubu, bsubv, bsubs, bsupv, bsupu, rzl_array, gc_ar
bsubumnc(:,1) = 0
bsubvmnc(:,1) = 0

! NOTE: This assumes that bsubs (from which bsubsmns is computed) is on the full grid.
! However, since crmn_o is passed as bsubs in fileout(),
! HERE, bsubs is ACTUALLY on the HALF-grid !!! WTF ???
! (The full-grid bsubs array is a LOCAL variable in jxbforce().)
! bsubs is on the full grid here: jxbforce puts it there.
bsubsmns(:,1) = 2*bsubsmns(:,2) - bsubsmns(:,3) ! extrapolation on full grid

bsupumnc(:,1) = 0
Expand Down
Loading