From c24881a08e17299b2b30bdbcbd865b7c8f957ef4 Mon Sep 17 00:00:00 2001 From: CharlesCNorton Date: Thu, 3 Sep 2026 04:29:55 -0400 Subject: [PATCH 1/2] Carry B_s to the wout file on the full radial mesh --- src/jxbforce.f90 | 9 +++++++-- src/read_wout_mod.f | 14 ++++++-------- src/wrout.f90 | 5 +---- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/jxbforce.f90 b/src/jxbforce.f90 index 6f3345da3..cce18d3a3 100644 --- a/src/jxbforce.f90 +++ b/src/jxbforce.f90 @@ -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$ @@ -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 @@ -521,6 +522,10 @@ SUBROUTINE jxbforce(bsupu, bsupv, bsubu, bsubv, bsubsh, & !bsubs(ns,:) = 2*bsubs(ns,:) - bsubs(ns-1,:) ! TODO: from ns, ns-1 to ns ??? bsubs(ns,:) = 2*bsubs(ns-1,:) - bsubs(ns-2,:) + ! 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 + diff --git a/src/read_wout_mod.f b/src/read_wout_mod.f index 7f647f02f..3c1d8dca0 100644 --- a/src/read_wout_mod.f +++ b/src/read_wout_mod.f @@ -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) @@ -744,8 +744,7 @@ 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) @@ -753,7 +752,7 @@ SUBROUTINE Compute_Currents(bsubsmnc_i, bsubsmns_i, 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 @@ -783,8 +782,7 @@ 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) @@ -792,7 +790,7 @@ SUBROUTINE Compute_Currents(bsubsmnc_i, bsubsmns_i, 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 diff --git a/src/wrout.f90 b/src/wrout.f90 index 67c169f67..7b36dfdd8 100644 --- a/src/wrout.f90 +++ b/src/wrout.f90 @@ -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 From 5d2cc78c6ea541ed62b65ab36ebefb7df5042d48 Mon Sep 17 00:00:00 2001 From: CharlesCNorton Date: Thu, 3 Sep 2026 05:11:47 -0400 Subject: [PATCH 2/2] Extrapolate B_s to the boundary from the outermost half point --- src/jxbforce.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/jxbforce.f90 b/src/jxbforce.f90 index cce18d3a3..6f49169b8 100644 --- a/src/jxbforce.f90 +++ b/src/jxbforce.f90 @@ -519,8 +519,9 @@ 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.