From a4a39571ce1bbcaab5eed7b4678de550c3b4759b Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Fri, 25 Oct 2024 10:50:58 +0200 Subject: [PATCH] src/pec.f90: don't rely on P_inside, I_cd_inside MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit When combining multiple independent simulations with `gray -s`, the statistics in sum-summary.txt are computed by summing the profiles pointwise and calling `postproc_profiles` with their (supposedly common) MHD equilibrium. The subroutine assumes the normalisations of dP/dV⋅dV and J_φ⋅dA are P_inside(ρ=1), I_cd_inside(ρ=1), however if the individual profiles have been computed with different versions of Gray, built with different compilers, or on different machines, the values may be very slightly off. Further, due to the use of the relation Δρ² = ⟨ρ²⟩ - ⟨ρ⟩², the profile widths are very sensible to the normalisation, producing wildly incorrect values. For example, due to the refactor of `magsurf_data`, a change of about 0.05% in dV resulted in Δρ² being off by 300%. --- src/gray_core.f90 | 4 ++-- src/main.f90 | 10 ++++---- src/pec.f90 | 60 +++++++++++++++++++++++++---------------------- 3 files changed, 39 insertions(+), 35 deletions(-) diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 98a6e91..e172a71 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -614,8 +614,8 @@ contains results%tables%ec_profiles, rhop_tab, rhot_tab, jphi_beam, & jcd_beam, dpdv_beam, currins_beam, pins_beam, index_rt) - call postproc_profiles(equil, pabs_beam, icd_beam, rhot_tab, dpdv_beam, & - jphi_beam, rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip, & + call postproc_profiles(equil, rhot_tab, dpdv_beam, & + jphi_beam, rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip, & rhotp, drhotp, rhotj, drhotj, dpdvmx, jphimx, ratjamx, ratjbmx) ! *compute profiles width for current beam if (results%tables%summary%active) then diff --git a/src/main.f90 b/src/main.f90 index 8c9a484..ef34ad6 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -279,8 +279,8 @@ contains end if ! Open combined output files - open(newunit=prof, file=output_dir//'/'//'sum-ec-profiles.txt', action='write') - open(newunit=summ, file=output_dir//'/'//'sum-summary.txt', action='write') + open(newunit=prof, file=output_dir//'/sum-ec-profiles.txt', action='write') + open(newunit=summ, file=output_dir//'/sum-summary.txt', action='write') ! Move to the filelist directory ! (the filepaths are assumed relative to it) @@ -358,9 +358,9 @@ contains ! Recompute the summary values call postproc_profiles(equil, & - results%Pabs, results%Icd, rhot_tab, results%dPdV, jphi, & - rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip, & - rhotp, drhotp, rhotj, drhotj, dpdvmx, jphimx, ratjamx, ratjbmx) + rhot_tab, results%dPdV, jphi, rhotpav, drhotpav, rhotjava, & + drhotjava, dpdvp, jphip, rhotp, drhotp, rhotj, drhotj, & + dpdvmx, jphimx, ratjamx, ratjbmx) write(fmtstr, '("(",i0,"(1x,e12.5),i5,7(1x,e12.5))")') nextracols+12 write(summ, trim(fmtstr)) & diff --git a/src/pec.f90 b/src/pec.f90 index 6b035b0..b191826 100644 --- a/src/pec.f90 +++ b/src/pec.f90 @@ -244,8 +244,8 @@ contains end subroutine pec_tab - subroutine postproc_profiles(equil, pabs, currt, rhot_tab, dpdv, ajphiv, & - rhotpav,drhotpav,rhotjava,drhotjava,dpdvp, ajphip, & + subroutine postproc_profiles(equil, rhot_tab, dpdv, ajphiv, & + rhotpav,drhotpav,rhotjava,drhotjava,dpdvp, ajphip, & rhotp, drhotp, rhotjfi, drhotjfi, dpdvmx,ajmxfi, ratjamx, ratjbmx) ! Radial average values over power and current density profile @@ -253,7 +253,6 @@ contains use gray_equil, only : abstract_equil class(abstract_equil), intent(in) :: equil - real(wp_), intent(in) :: pabs,currt real(wp_), intent(in) :: rhot_tab(:) real(wp_), intent(in) :: dpdv(:), ajphiv(:) real(wp_), intent(out) :: rhotpav, drhotpav, dpdvp @@ -262,39 +261,44 @@ contains real(wp_), intent(out) :: rhotjfi, drhotjfi, ajmxfi real(wp_), intent(out) :: ratjamx, ratjbmx - real(wp_) :: sccsa, ratjplmx, rhopjava, rhoppav - real(wp_) :: rhotjav, rhot2pav, rhot2java, dvdrhotav, dadrhotava + real(wp_) :: P_tot, I_tot, ratjplmx, rhopjava, rhoppav + real(wp_) :: rhotjav, dvdrhotav, dadrhotava - rhotpav = 0 - rhot2pav = 0 - rhotjav = 0 - rhotjava = 0 - rhot2java = 0 - - if (pabs > 0) then - rhotpav = sum(rhot_tab *dpdv*dvol)/pabs - rhot2pav = sum(rhot_tab**2*dpdv*dvol)/pabs + ! ⟨ρ_t⟩, Δρ_t using dP/dV⋅dV as a PDF + ! Note: the factor √8 is to match the FWHM of a Gaussian + P_tot = sum(dpdv*dvol) + if (P_tot > 0) then + rhotpav = sum(rhot_tab * dpdv*dvol)/P_tot + drhotpav = sqrt(8 * sum((rhot_tab - rhotpav)**2 * dpdv*dvol)/P_tot) + else + rhotpav = 0 + drhotpav = 0 end if - if (abs(currt) > 0) then - rhotjav = sum(rhot_tab*ajphiv*darea)/currt - end if - sccsa = sum(abs(ajphiv)*darea) - if (sccsa > 0) then - rhotjava = sum(rhot_tab *abs(ajphiv)*darea)/sccsa - rhot2java = sum(rhot_tab**2*abs(ajphiv)*darea)/sccsa + ! ⟨ρ_t⟩ using J_φ⋅dA as a PDF + I_tot = sum(ajphiv*darea) + if (abs(I_tot) > 0) then + rhotjav = sum(rhot_tab*ajphiv*darea)/I_tot + else + rhotjav = 0 end if - ! factor sqrt(8) = 2sqrt(2) to match full width of gaussian profile - drhotpav = sqrt(8._wp_*(max(rhot2pav -rhotpav**2,comp_eps))) - drhotjava = sqrt(8._wp_*(max(rhot2java-rhotjava**2,comp_eps))) + ! ⟨ρ_t⟩, Δρ_t using |J_φ|⋅dA as a PDF + I_tot = sum(abs(ajphiv)*darea) + if (I_tot > 0) then + rhotjava = sum(rhot_tab * abs(ajphiv)*darea)/I_tot + drhotjava = sqrt(8* sum((rhot_tab - rhotjava)**2 * abs(ajphiv)*darea)/I_tot) + else + rhotjava = 0 + drhotjava = 0 + end if rhoppav = equil%tor2pol(rhotpav) rhopjava = equil%tor2pol(rhotjava) - if (pabs > 0) then + if (P_tot > 0) then call equil%flux_average(rhoppav, dVdrho_t=dvdrhotav) - dpdvp = pabs*2.0_wp_/(sqrt(pi)*drhotpav*dvdrhotav) + dpdvp = P_tot * 2/(sqrt(pi)*drhotpav*dvdrhotav) call profwidth(rhot_tab,dpdv,rhotp,dpdvmx,drhotp) else dpdvp = 0 @@ -303,10 +307,10 @@ contains drhotp = 0 end if - if (sccsa > 0) then + if (I_tot > 0) then call equil%flux_average(rhopjava, dAdrho_t=dadrhotava, ratio_astra_tor=ratjamx, & ratio_jintrac_tor=ratjbmx, ratio_paral_tor=ratjplmx) - ajphip = currt*2.0_wp_/(sqrt(pi)*drhotjava*dadrhotava) + ajphip = I_tot * 2/(sqrt(pi)*drhotjava*dadrhotava) call profwidth(rhot_tab,ajphiv,rhotjfi,ajmxfi,drhotjfi) else ajphip = 0