Merge pull request #23 from ProjectTorreyPines/handle_gsi_change

Handle gsi change when extended grid displaced B2.5 grid

Merging as ProjectTorreyPines/SOLPS2imas.jl#30 as been merged.

samples/time_dep_edge_profiles_with_equilibrium.json.dvc

@@ -1,12 +1,12 @@
-md5: c323be561df2fd7057b611798041547c
+md5: 8d18819e21169fdfbe55c13be90a693a
 frozen: true
 deps:
 - path: ITER_Lore_2296_00000/IMAS/time_dep_edge_profiles_with_equilibrium.json
   repo:
     url: git@github.com:ProjectTorreyPines/SOLPSTestSamples.git
-    rev_lock: 1593ec15e791e8d1987a3aa6bafbf0a8968e2f08
+    rev_lock: 6b5d423c56e49da55a4822ad07d0f525c5c53038
 outs:
-- md5: a1848b5908bf756323fb3d03a2d2632c
-  size: 40580800
+- md5: 8fbd26ff4b41c8f56a63f0918bb554dc
+  size: 40872278
   hash: md5
   path: time_dep_edge_profiles_with_equilibrium.json

src/interferometer.jl

@@ -43,15 +43,20 @@ line integrated electron density if not present
 function add_interferometer!(
     config::String=default_ifo,
     @nospecialize(ids::IMASDD.dd)=IMASDD.dd();
-    overwrite::Bool=false, verbose::Bool=false,
+    overwrite::Bool=false, verbose::Bool=false, rtol::Float64=1e-3, n_e_gsi::Int=5,
 )::IMASDD.dd
     if endswith(config, ".json")
         config_dict = convert_strings_to_symbols(IMASDD.JSON.parsefile(config))
-        add_interferometer!(config_dict, ids; overwrite, verbose)
+        add_interferometer!(
+            config_dict,
+            ids;
+            overwrite=overwrite,
+            verbose=verbose,
+            n_e_gsi=n_e_gsi,
+        )
     else
         error("Only JSON files are supported.")
     end
-    compute_interferometer(ids)
     return ids
 end
 
@@ -68,7 +73,7 @@ electron density if not present
 function add_interferometer!(
     config::Dict{Symbol, Any},
     @nospecialize(ids::IMASDD.dd)=IMASDD.dd();
-    overwrite::Bool=false, verbose::Bool=false,
+    overwrite::Bool=false, verbose::Bool=false, rtol::Float64=1e-3, n_e_gsi::Int=5,
 )::IMASDD.dd
     # Check for duplicates
     if length(ids.interferometer.channel) > 0
@@ -113,7 +118,7 @@ function add_interferometer!(
             )
     end
     IMASDD.dict2imas(config, ids; verbose)
-    compute_interferometer(ids)
+    compute_interferometer(ids; rtol=rtol, n_e_gsi=n_e_gsi)
     return ids
 end
 
@@ -123,29 +128,46 @@ end
   - Calculate phase_to_n_e_line if not present for each wavelength
   - Compute the line integrated electron density if not present
 """
-function compute_interferometer(@nospecialize(ids::IMASDD.dd), rtol::Float64=1e-3)
+function compute_interferometer(
+    @nospecialize(ids::IMASDD.dd);
+    rtol::Float64=1e-3,
+    n_e_gsi::Int=5,
+)
+    # Compute phase_to_n_e_line if not present
+    for ch ∈ ids.interferometer.channel
+        k = @SVector[2π / ch.wavelength[ii].value for ii ∈ 1:2]
+        for i1 ∈ 1:2
+            lam = ch.wavelength[i1]
+            i2 = i1 % 2 + 1
+            if IMASDD.ismissing(lam, :phase_to_n_e_line)
+                # Taken from https://doi.org/10.1063/1.1138037
+                lam.phase_to_n_e_line =
+                    (
+                        2 * m_e * ε_0 * c_0^2 / e^2 * (k[i1] * k[i2]^2) /
+                        (k[i2]^2 - k[i1]^2)
+                    ).val
+            end
+        end
+    end
+
     fix_eq_time_idx = length(ids.equilibrium.time_slice) == 1
     fix_ep_grid_ggd_idx = length(ids.edge_profiles.grid_ggd) == 1
 
     ep_grid_ggd = ids.edge_profiles.grid_ggd[1]
     ep_space = ep_grid_ggd.space[1]
     sep_bnd = get_sep_bnd(ep_grid_ggd)
-    # Using -5 for now to use the SOLPS edge profile grid only
-    TPS_mats = get_TPS_mats(ep_grid_ggd, -5)
 
     epggd = ids.edge_profiles.ggd
     cpp1d = ids.core_profiles.profiles_1d
     nt = length(epggd)
 
+    TPS_mats = get_TPS_mats(ep_grid_ggd, n_e_gsi)
+
     ep_n_e_list = [
         interp(
             epggd[ii].electrons.density,
-            update_TPS_mats(ii, fix_ep_grid_ggd_idx, ids, -5, TPS_mats),
-            5,
-            # Note the grid_subset index 5 is used here to get the electron density
-            # data from SOLPS edge profile. This is a bug in SD4SOLPS. On adding
-            # edge extension, it should update all quantitites that refered to
-            # grid_subset index 5 to -5.
+            update_TPS_mats(ii, fix_ep_grid_ggd_idx, ids, n_e_gsi, TPS_mats),
+            n_e_gsi,
         ) for ii ∈ eachindex(epggd)
     ]
     cp_n_e_list = [
@@ -157,41 +179,35 @@ function compute_interferometer(@nospecialize(ids::IMASDD.dd), rtol::Float64=1e-
     ]
 
     for ch ∈ ids.interferometer.channel
-        k = @SVector[2π / ch.wavelength[ii].value for ii ∈ 1:2]
-        for i1 ∈ 1:2
-            lam = ch.wavelength[i1]
-            i2 = i1 % 2 + 1
-            if IMASDD.ismissing(lam, :phase_to_n_e_line)
-                # Taken from https://doi.org/10.1063/1.1138037
-                lam.phase_to_n_e_line =
-                    (
-                        2 * m_e * ε_0 * c_0^2 / e^2 * (k[i1] * k[i2]^2) /
-                        (k[i2]^2 - k[i1]^2)
-                    ).val
-            end
+        # Special case when measurement has not been made for some time steps 
+        # but edge profile data exists
+        proceed =
+            (IMASDD.ismissing(ch.n_e_line, :time) || IMASDD.isempty(ch.n_e_line.time))
+        if !proceed
+            proceed = length(ch.n_e_line.time) < length(epggd)
         end
-        # Special case when measurement has not been made but edge profile data exists
-        if (IMASDD.ismissing(ch.n_e_line, :time) || IMASDD.isempty(ch.n_e_line.time)) &&
-           length(ids.edge_profiles.ggd) > 0
-            ch.n_e_line.time = zeros(nt)
-            ch.n_e_line_average.time = zeros(nt)
-            ch.n_e_line.data = zeros(nt)
-            ch.n_e_line_average.data = zeros(nt)
-            for lam ∈ ch.wavelength
-                lam.phase_corrected.time = zeros(nt)
-                lam.phase_corrected.data = zeros(nt)
-            end
-
+        if proceed
             # Check if core_profile is available
             if length(cpp1d) != length(epggd)
                 error(
-                    "Number of edge profiles time slices does not match number of " \
-                    "core profile time slices. Please ensure core profile data for " \
-                    "electron density is present in the data structure. You might " \
-                    "want to run SD4SOLPS.fill_in_extrapolated_core_profile!" \
+                    "Number of edge profiles time slices does not match number " \
+                    "of core profile time slices. Please ensure core profile " \
+                    "data for electron density is present in the data structure. " \
+                    "You might want to run " \
+                    "SD4SOLPS.fill_in_extrapolated_core_profile!" \
                     "(dd, \"electrons.density\"; method=core_method))",
                 )
             end
+
+            resize!(ch.n_e_line.time, nt)
+            resize!(ch.n_e_line_average.time, nt)
+            resize!(ch.n_e_line.data, nt)
+            resize!(ch.n_e_line_average.data, nt)
+            for lam ∈ ch.wavelength
+                resize!(lam.phase_corrected.time, nt)
+                resize!(lam.phase_corrected.data, nt)
+            end
+
             # Parametrize line of sight
             fp = rzphi2xyz(ch.line_of_sight.first_point)
             sp = rzphi2xyz(ch.line_of_sight.second_point)
@@ -202,7 +218,8 @@ function compute_interferometer(@nospecialize(ids::IMASDD.dd), rtol::Float64=1e-
             else
                 chord_points = (fp, sp, tp)
             end
-            core_chord_length = get_core_chord_length(sep_bnd, ep_space, chord_points)
+            core_chord_length =
+                get_core_chord_length(sep_bnd, ep_space, chord_points)
 
             for ii ∈ eachindex(epggd)
                 ch.n_e_line.time[ii] = epggd[ii].time
@@ -230,8 +247,9 @@ function compute_interferometer(@nospecialize(ids::IMASDD.dd), rtol::Float64=1e-
                         )
                     end
 
-                ch.n_e_line.data[ii] = quadgk(integ, 0, 1; rtol)[1]
-                ch.n_e_line_average.data[ii] = ch.n_e_line.data[ii] / core_chord_length
+                ch.n_e_line.data[ii] = quadgk(integ, 0, 1; rtol=rtol)[1]
+                ch.n_e_line_average.data[ii] =
+                    ch.n_e_line.data[ii] / core_chord_length
                 for lam ∈ ch.wavelength
                     lam.phase_corrected.time[ii] = epggd[ii].time
                     lam.phase_corrected.data[ii] =

src/langmuir_probes.jl

@@ -28,16 +28,16 @@ function compute_langmuir_probes(
     ids::IMASDD.dd;
     ne_noise::Union{Noise, Nothing}=nothing,
     te_noise::Union{Noise, Nothing}=nothing,
+    n_e_gsi::Int=5,
 )
     epggd = ids.edge_profiles.ggd
     nt = length(epggd)
     fix_ep_grid_ggd_idx = length(ids.edge_profiles.grid_ggd) == 1
     ep_grid_ggd = ids.edge_profiles.grid_ggd[1]
-    # Using -5 for now to use the SOLPS edge profile grid only
-    TPS_mats_all_cells = get_TPS_mats(ep_grid_ggd, -5)
+    TPS_mats_all_cells = get_TPS_mats(ep_grid_ggd, n_e_gsi)
     # Get the edge profile interpolation functions
-    ep_n_e = interp(epggd[1].electrons.density, TPS_mats_all_cells, 5)
-    ep_t_e = interp(epggd[1].electrons.temperature, TPS_mats_all_cells, 5)
+    ep_n_e = interp(epggd[1].electrons.density, TPS_mats_all_cells, n_e_gsi)
+    ep_t_e = interp(epggd[1].electrons.temperature, TPS_mats_all_cells, n_e_gsi)
 
     # Initialize langmuir probe data
     init_data!.(ids.langmuir_probes.embedded, nt)
@@ -46,10 +46,11 @@ function compute_langmuir_probes(
         if !fix_ep_grid_ggd_idx && ii > 1
             # If grid_ggd is evolving with time, update boundaries
             ep_grid_ggd = ids.edge_profiles.grid_ggd[ii]
-            TPS_mats_all_cells = get_TPS_mats(ep_grid_ggd, -5)
+            TPS_mats_all_cells = get_TPS_mats(ep_grid_ggd, n_e_gsi)
             # Update the edge profile interpolation functions
-            ep_n_e = interp(epggd[ii].electrons.density, TPS_mats_all_cells, 5)
-            ep_t_e = interp(epggd[ii].electrons.temperature, TPS_mats_all_cells, 5)
+            ep_n_e = interp(epggd[ii].electrons.density, TPS_mats_all_cells, n_e_gsi)
+            ep_t_e =
+                interp(epggd[ii].electrons.temperature, TPS_mats_all_cells, n_e_gsi)
         end
 
         for emb_lp ∈ ids.langmuir_probes.embedded

test/runtests.jl

@@ -31,7 +31,11 @@ if args["interferometer"]
             json2imas(
                 "$(@__DIR__)/../samples/time_dep_edge_profiles_with_equilibrium.json",
             )
-        add_interferometer!("$(@__DIR__)/../src/default_interferometer.json", ids)
+        add_interferometer!(
+            "$(@__DIR__)/../src/default_interferometer.json",
+            ids;
+            n_e_gsi=-5,
+        )
         # Just checking if the function runs through for now
         for ch ∈ ids.interferometer.channel
             println()
@@ -56,7 +60,7 @@ if args["interferometer"]
         try
             add_interferometer!(
                 "$(@__DIR__)/../samples/test_interferometer_same_names.json",
-                ids,
+                ids; n_e_gsi=-5,
             )
         catch err
             @test isa(err, OverwriteAttemptError)
@@ -65,7 +69,7 @@ if args["interferometer"]
         add_interferometer!(
             "$(@__DIR__)/../samples/test_interferometer_same_names.json",
             ids;
-            overwrite=true,
+            overwrite=true, n_e_gsi=-5,
         )
         @test ids.interferometer.channel[1].name == "V1"
         @test ids.interferometer.channel[1].line_of_sight.first_point.r == 5.2
@@ -74,6 +78,7 @@ if args["interferometer"]
         # Try appending new interfermeter channels
         add_interferometer!(
             "$(@__DIR__)/../samples/test_interferometer_new_channels.json", ids;
+            n_e_gsi=-5,
         )
         for ch ∈ ids.interferometer.channel
             println()
@@ -114,7 +119,8 @@ if args["langmuir_probes"]
         add_langmuir_probes!(
             "$(@__DIR__)/../src/default_langmuir_probes.json",
             ids;
-            ne_noise,
+            ne_noise=ne_noise,
+            n_e_gsi=-5,
         )
         # Just checking if the function runs through for now
         for lp ∈ ids.langmuir_probes.embedded