diff --git a/samples/time_dep_edge_profiles_with_equilibrium.json.dvc b/samples/time_dep_edge_profiles_with_equilibrium.json.dvc index 181bd04..d33a57c 100644 --- a/samples/time_dep_edge_profiles_with_equilibrium.json.dvc +++ b/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 diff --git a/src/interferometer.jl b/src/interferometer.jl index 092a2c2..df49039 100644 --- a/src/interferometer.jl +++ b/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] = diff --git a/src/langmuir_probes.jl b/src/langmuir_probes.jl index 26af3ca..0ebdde5 100644 --- a/src/langmuir_probes.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index ccb8284..458681c 100644 --- a/test/runtests.jl +++ b/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