Update fluxsurfaces.jl

src/physics/fast.jl

@@ -180,7 +180,7 @@ function critical_energy(cp1d::IMAS.core_profiles__profiles_1d, mf::Real, Zf::Re
     Te = cp1d.electrons.temperature
 
     avg_cmr = sum(ion.density_thermal .* (avgZ(ion.element[1].z_n, sum(ion.temperature) / length(ion.temperature))^2) / ion.element[1].a for ion in cp1d.ion) ./ ne
-    Ec = 14.8 .* mf .* Te .* avg_cmr .^ 1.5
+    Ec = 14.8 .* mf .* Te .* avg_cmr .^ (2.0 / 3.0)
     if !approximate
         Ec1 = critical_energy(cp1d, 1, mf, Zf; approximate)
         Ec = Ec .* (Ec1 ./ Ec[1])

src/physics/fluxsurfaces.jl

@@ -1765,7 +1765,12 @@ function flux_surfaces(eqt::equilibrium__time_slice{T1}, wall_r::AbstractVector{
 
     # volume
     eqt1d.volume = cumtrapz(eqt1d.psi, eqt1d.dvolume_dpsi)
-    eqt.global_quantities.volume = eqt1d.volume[end]
+    # int_fluxexpansion_dl is inaccurate at the LCFS near X-points (Bp→0 diverges numerically),
+    # so correct dvolume_dpsi[end] and volume[end] using the geometrically exact boundary outline volume
+    vol_lcfs = π * abs(trapz(eqt.boundary.outline.z, eqt.boundary.outline.r .^ 2))
+    eqt1d.dvolume_dpsi[end] = 2 * (vol_lcfs - eqt1d.volume[end-1]) / (eqt1d.psi[end] - eqt1d.psi[end-1]) - eqt1d.dvolume_dpsi[end-1]
+    eqt1d.volume[end] = vol_lcfs
+    eqt.global_quantities.volume = vol_lcfs
 
     # phi
     tmp .= eqt1d.f .* eqt1d.gm1 / (2π)

src/physics/sources.jl

@@ -284,7 +284,7 @@ function total_power_inside(
 )
     exclude_indexes = Int[]
     if !include_time_derivative
-        exclude_indexes = [11]
+        exclude_indexes = [409]
     else
         exclude_indexes = Int[]
     end