fix handling of pedestal density

ProjectTorreyPines · Feb 12, 2025 · d7c6b6c · d7c6b6c
1 parent a5bd6f3
commit d7c6b6c
Show file tree

Hide file tree

Showing 3 changed files with 20 additions and 20 deletions.
diff --git a/src/actors/pedestal/pedestal_actor.jl b/src/actors/pedestal/pedestal_actor.jl
@@ -57,8 +57,10 @@ end
 function ActorPedestal(dd::IMAS.dd, par::FUSEparameters__ActorPedestal, act::ParametersAllActors; kw...)
     logging_actor_init(ActorPedestal)
     par = par(kw...)
-    eped_actor = ActorEPED(dd, act.ActorEPED; par.rho_nml, par.rho_ped, par.T_ratio_pedestal, par.Te_sep, par.ip_from, par.βn_from, ne_from=:core_profiles, zeff_from=:core_profiles)
-    wped_actor = ActorWPED(dd, act.ActorWPED; par.rho_nml, par.rho_ped, par.T_ratio_pedestal, par.Te_sep, par.ip_from, par.βn_from, ne_from=:core_profiles, zeff_from=:core_profiles)
+    eped_actor =
+        ActorEPED(dd, act.ActorEPED; par.rho_nml, par.rho_ped, par.T_ratio_pedestal, par.Te_sep, par.ip_from, par.βn_from, ne_from=:core_profiles, zeff_from=:core_profiles)
+    wped_actor =
+        ActorWPED(dd, act.ActorWPED; par.rho_nml, par.rho_ped, par.T_ratio_pedestal, par.Te_sep, par.ip_from, par.βn_from, ne_from=:core_profiles, zeff_from=:core_profiles)
     noop = ActorNoOperation(dd, act.ActorNoOperation)
     actor = ActorPedestal(dd, par, act, noop, wped_actor, eped_actor, noop, noop, Symbol[], -Inf, -Inf, -Inf, IMAS.core_profiles__profiles_1d())
     actor.replay_actor = ActorReplay(dd, act.ActorReplay, actor)
@@ -194,17 +196,19 @@ function pedestal_density_tanh(dd::IMAS.dd, par::FUSEparameters__ActorPedestal)
     rho = cp1d.grid.rho_tor_norm
 
     rho09 = 0.9
-    w_ped = 0.05
+    w_ped_ne = 0.05
     ne_ped_old = IMAS.interp1d(rho, cp1d.electrons.density_thermal).(rho09)
     ne_ped = IMAS.get_from(dd, Val{:ne_ped}, par.ne_from, rho09)
     cp1d.electrons.density_thermal[end] = ne_ped / 4.0
-    cp1d.electrons.density_thermal = IMAS.blend_core_edge_Hmode(cp1d.electrons.density_thermal, rho, ne_ped, w_ped, par.rho_nml, par.rho_ped)
+    ne = IMAS.blend_core_edge_Hmode(cp1d.electrons.density_thermal, rho, ne_ped, w_ped_ne, par.rho_nml, par.rho_ped)
+    cp1d.electrons.density_thermal = IMAS.ped_height_at_09(rho, ne, ne_ped)
     for ion in cp1d.ion
         if !ismissing(ion, :density_thermal)
             ni_ped_old = IMAS.interp1d(rho, ion.density_thermal).(rho09)
-            ni_ped = ni_ped_old  / ne_ped_old * ne_ped
+            ni_ped = ni_ped_old / ne_ped_old * ne_ped
             ion.density_thermal[end] = ni_ped / 4.0
-            ion.density_thermal = IMAS.blend_core_edge_Hmode(ion.density_thermal, rho, ni_ped, w_ped, par.rho_nml, par.rho_ped)
+            ni = IMAS.blend_core_edge_Hmode(ion.density_thermal, rho, ni_ped, w_ped_ne, par.rho_nml, par.rho_ped)
+            ion.density_thermal = IMAS.ped_height_at_09(rho, ni, ni_ped)
         end
     end
 end

diff --git a/src/actors/transport/eped_profiles_actor.jl b/src/actors/transport/eped_profiles_actor.jl
@@ -67,7 +67,7 @@ function _step(actor::ActorEPEDprofiles)
     tval = IMAS.Hmode_profiles(Te[end], Te_ped, Te[1], length(cp1d.grid.rho_tor_norm), par.Te_shaping, par.Te_shaping, wped)
     cp1d.electrons.temperature = tval
     if any(tval .< 0)
-        error("Te profile is negative for T0=$(Te[1]) eV and Tped=$(Te_ped) eV")
+        error("Te profile is negative for T0=$(Te[1]) [eV] and Tped=$(Te_ped) [eV]")
     end
 
     # update ion temperature profiles using:
@@ -82,7 +82,7 @@ function _step(actor::ActorEPEDprofiles)
         par.Te_shaping,
         wped)
     if any(tval_ions .< 0)
-        error("Ti profile is negative for T0=$(Te[1]*par.T_ratio_core) eV and Tped=$(Te_ped*par.T_ratio_pedestal) eV")
+        error("Ti profile is negative for T0=$(Te[1]*par.T_ratio_core) [eV] and Tped=$(Te_ped*par.T_ratio_pedestal) [eV]")
     end
     for ion in cp1d.ion
         ion.temperature = tval_ions
@@ -97,10 +97,10 @@ function _step(actor::ActorEPEDprofiles)
     for (ii, ion) in enumerate(cp1d.ion)
         ion_fractions[ii, :] = ion.density_thermal ./ ne
     end
-    nval = IMAS.Hmode_profiles(ne[end], ne_ped, ne[1], length(cp1d.grid.rho_tor_norm), par.ne_shaping, par.ne_shaping, wped)
-    cp1d.electrons.density_thermal = nval
+    nval = IMAS.Hmode_profiles(ne[end], ne_ped, ne[1], length(cp1d.grid.rho_tor_norm), par.ne_shaping, par.ne_shaping, 0.05)
+    cp1d.electrons.density_thermal = IMAS.ped_height_at_09(cp1d.grid.rho_tor_norm, nval, ne_ped)
     if any(nval .< 0)
-        error("ne profile is negative for n0=$(ne[1]) /m^3 and Tped=$(ne_ped) /m^3")
+        error("ne profile is negative for n0=$(ne[1]) [m⁻³] and ne_ped=$(ne_ped) [m⁻³]")
     end
 
     # update ion density profiles using

diff --git a/src/ddinit/init_core_profiles.jl b/src/ddinit/init_core_profiles.jl
@@ -133,21 +133,17 @@ function init_core_profiles!(
     @assert ne_core_to_ped_ratio > 1.0
     @assert plasma_mode in (:H_mode, :L_mode)
 
-    if plasma_mode == :L_mode
-        ne_core_to_ped_ratio = 1.0 + 1.0 / w_ped * (1.0 - w_ped)
-    end
-
     cp1d = resize!(cp.profiles_1d)
 
     # grid
     cp1d.grid.rho_tor_norm = range(0, 1, ngrid)
 
     # Density
-    if plasma_mode == :H_mode
-        cp1d.electrons.density_thermal = IMAS.Hmode_profiles(ne_sep_to_ped_ratio, 1.0, ne_core_to_ped_ratio, ngrid, ne_shaping, ne_shaping, w_ped)
-    else
-        cp1d.electrons.density_thermal = IMAS.Lmode_profiles(ne_sep_to_ped_ratio, 1.0, ne_core_to_ped_ratio, ngrid, ne_shaping, 1.0, w_ped)
-    end
+    # first we start with a "unit" density profile...
+    w_ped_ne = 0.05
+    cp1d.electrons.density_thermal = IMAS.Hmode_profiles(ne_sep_to_ped_ratio, 1.0, ne_core_to_ped_ratio, ngrid, ne_shaping, ne_shaping, w_ped_ne)
+    cp1d.electrons.density_thermal = IMAS.ped_height_at_09(cp1d.grid.rho_tor_norm, cp1d.electrons.density_thermal, 1.0)
+    # ...which then we scale according to :ne_setting and :ne_value
     if ne_setting == :ne_ped
         ne_ped = ne_value
     elseif ne_setting == :greenwald_fraction_ped