From 37d4f304072271a1fdf8c01002610ef05b0c9266 Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Tue, 26 Sep 2023 13:31:48 -0700 Subject: [PATCH 01/12] Add formatting --- .JuliaFormatter.toml | 27 ++++++++++++++++++++++ .github/workflows/format_check.yml | 36 ++++++++++++++++++++++++++++++ README.md | 7 +++++- 3 files changed, 69 insertions(+), 1 deletion(-) create mode 100644 .JuliaFormatter.toml create mode 100644 .github/workflows/format_check.yml diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 0000000..aa7f0eb --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1,27 @@ +align_assignment = true +align_conditional = true +align_matrix = true +align_pair_arrow = true +align_struct_field = true +always_for_in = true +always_use_return = true +annotate_untyped_fields_with_any = false +conditional_to_if = true +for_in_replacement = "∈" +format_docstrings = true +import_to_using = true +indent = 4 +indent_submodule = true +join_lines_based_on_source = true +long_to_short_function_def = true +margin = 88 +normalize_line_endings = "unix" +pipe_to_function_call = true +remove_extra_newlines = true +separate_kwargs_with_semicolon = true +surround_whereop_typeparameters = true +trailing_comma = true +whitespace_in_kwargs = false +whitespace_ops_in_indices = false +whitespace_typedefs = true +yas_style_nesting = true diff --git a/.github/workflows/format_check.yml b/.github/workflows/format_check.yml new file mode 100644 index 0000000..4e46a44 --- /dev/null +++ b/.github/workflows/format_check.yml @@ -0,0 +1,36 @@ +name: Format Check + +on: + push: + branches: ["master", "dev", "format"] + pull_request: + branches: ["master", "dev"] +jobs: + check: + runs-on: ${{ matrix.os }} + strategy: + matrix: + julia-version: [1.9.3] + julia-arch: [x86] + os: [ubuntu-latest] + steps: + - uses: julia-actions/setup-julia@latest + with: + version: ${{ matrix.julia-version }} + + - uses: actions/checkout@v1 + - name: Install JuliaFormatter and format + run: | + julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))' + julia -e 'using JuliaFormatter; format(".", verbose=true)' + - name: Format check + run: | + julia -e ' + out = Cmd(`git diff --name-only`) |> read |> String + if out == "" + exit(0) + else + @error "Some files have not been formatted !!!" + write(stdout, out) + exit(1) + end' diff --git a/README.md b/README.md index 09444f2..dafdc7a 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,7 @@ # SynthDiag.jl -Package that defined synthetic diagnostics to generate sensor data based on plasma profile. + +![Format Check](https://github.com/ProjectTorreyPines/GGDUtils.jl/actions/workflows/format_check.yml/badge.svg) + + +Package that defined synthetic diagnostics to generate sensor data based on plasma +profile. From b3a3d3f54fcd4c650f9628e2c2eb405e40424847 Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Tue, 26 Sep 2023 14:02:49 -0700 Subject: [PATCH 02/12] Initializing package --- .gitignore | 5 +++++ Project.toml | 4 ++++ src/SynthDiag.jl | 3 +++ 3 files changed, 12 insertions(+) create mode 100644 .gitignore create mode 100644 Project.toml create mode 100644 src/SynthDiag.jl diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..6b72f03 --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +# File generated by Pkg, the package manager, based on a corresponding Project.toml +# It records a fixed state of all packages used by the project. As such, it should not be +# committed for packages, but should be committed for applications that require a static +# environment. +Manifest.toml diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..35acc32 --- /dev/null +++ b/Project.toml @@ -0,0 +1,4 @@ +name = "SynthDiag" +uuid = "c6e34158-aa22-49e4-bb12-15cbcb518ce6" +authors = ["Anchal Gupta "] +version = "0.1.0" diff --git a/src/SynthDiag.jl b/src/SynthDiag.jl new file mode 100644 index 0000000..3bbcee9 --- /dev/null +++ b/src/SynthDiag.jl @@ -0,0 +1,3 @@ +module SynthDiag + +end # module SynthDiag From 9324622ed8f23ca81eed783f8cfd25b5c8ea60d6 Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Tue, 26 Sep 2023 18:03:02 -0700 Subject: [PATCH 03/12] WIP Adding reading of interferometer configuration and computation roadmap --- src/SynthDiag.jl | 4 ++ src/default_interferometer.yml | 82 +++++++++++++++++++++++++++++ src/interferometer.jl | 94 ++++++++++++++++++++++++++++++++++ 3 files changed, 180 insertions(+) create mode 100644 src/default_interferometer.yml create mode 100644 src/interferometer.jl diff --git a/src/SynthDiag.jl b/src/SynthDiag.jl index 3bbcee9..fbd4a4b 100644 --- a/src/SynthDiag.jl +++ b/src/SynthDiag.jl @@ -1,3 +1,7 @@ module SynthDiag +using OMAS: OMAS + +include("$(@__DIR__)/interferometer.jl") + end # module SynthDiag diff --git a/src/default_interferometer.yml b/src/default_interferometer.yml new file mode 100644 index 0000000..55c0d50 --- /dev/null +++ b/src/default_interferometer.yml @@ -0,0 +1,82 @@ +interferometer: + channel: + - name: V1 + identifier: 0 + line_of_sight: + first_point: + phi: 0.0 + r: 5.0 + z: -5.0 + second_point: + phi: 0.0 + r: 5.0 + z: 5.0 + third_point: + phi: 0.0 + r: 5.0 + z: -5.0 + wavelength: + - phase_to_n_e_line: + value: 10.6e-6 # in m, for CO2 + - phase_to_n_e_line: + value: 0.633e-6 # in m, for HeNe + - name: V2 + identifier: 1 + line_of_sight: + first_point: + phi: 0.0 + r: 6.0 + z: -5.0 + second_point: + phi: 0.0 + r: 6.0 + z: 5.0 + third_point: + phi: 0.0 + r: 6.0 + z: -5.0 + wavelength: + - phase_to_n_e_line: + value: 10.6e-6 # in m, for CO2 + - phase_to_n_e_line: + value: 0.633e-6 # in m, for HeNe + - name: V3 + identifier: 1 + line_of_sight: + first_point: + phi: 0.0 + r: 7.0 + z: -5.0 + second_point: + phi: 0.0 + r: 7.0 + z: 5.0 + third_point: + phi: 0.0 + r: 7.0 + z: -5.0 + wavelength: + - phase_to_n_e_line: + value: 10.6e-6 # in m, for CO2 + - phase_to_n_e_line: + value: 0.633e-6 # in m, for HeNe + - name: H1 + identifier: 1 + line_of_sight: + first_point: + phi: 0.0 + r: 9.0 + z: 0.0 + second_point: + phi: 0.0 + r: 3.0 + z: 0.0 + third_point: + phi: 0.0 + r: 9.0 + z: 0.0 + wavelength: + - phase_to_n_e_line: + value: 10.6e-6 # in m, for CO2 + - phase_to_n_e_line: + value: 0.633e-6 # in m, for HeNe diff --git a/src/interferometer.jl b/src/interferometer.jl new file mode 100644 index 0000000..bb00dc2 --- /dev/null +++ b/src/interferometer.jl @@ -0,0 +1,94 @@ +import PhysicalConstants.CODATA2018: c_0, ε_0, m_e, e +import SD4SOLPS: fill_in_extrapolated_core_profile, check_rho_1d, add_rho_to_equilibrium + +default_ifo = "$(@__DIR__)/default_interferometer.yml" + +""" + add_interferometer( + @nospecialize(ids::OMAS.dd)=OMAS.dd(), + config::String=default_ifo, + +)::OMAS.dd + +Add interferometer to IMAS structure using either a YAML or JSON file and compute the +line integrated electron density if not present +""" +function add_interferometer( + @nospecialize(ids::OMAS.dd)=OMAS.dd(), + config::String=default_ifo, +)::OMAS.dd + if endswith(config, ".yml") || endswith(config, ".yaml") + OMAS.yaml2imas(config, ids) + elseif endswith(config, ".json") + OMAS.json2imas(config, ids) + end + compute_interferometer(ids) + return ids +end + +""" + add_interferometer( + @nospecialize(ids::OMAS.dd)=OMAS.dd(), + config::Dict{Symbol, Any}, + +)::OMAS.dd + +Add interferometer to IMAS structure using a Dict and compute the line integrated +electron density if not present +""" +function add_interferometer( + @nospecialize(ids::OMAS.dd)=OMAS.dd(), + config::Dict{Symbol, Any}, +)::OMAS.dd + OMAS.dict2imas(config, ids) + compute_interferometer(ids) + return ids +end + +""" + compute_interferometer(@nospecialize(ids::OMAS.dd)) + + - 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::OMAS.dd)) + for ch ∈ ids.interferometer.channel + k = [2 * π / ch.wavelength[ii].value for ii ∈ 1:2] + for i1 ∈ 1:2 + lam = ch.wavelength[i1] + i2 = i1 % 2 + 1 + if lam.phase_to_n_e_line == 0 + # 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) + end + end + # Special case when measurement has not been made but edge profile data exists + if length(ch.n_e_line.time) == 0 && length(ids.edge_profiles.ggd) > 0 + epggd = ids.edge_profiles.ggd + cpp1d = ids.core_profiles.profiles_1d + ch.n_e_line.time = ch.n_e_line.value = zeros(length(epggd)) + # Check if equilibrium data contains rho, if not add it + if !check_rho_1d(ids; time_slice=1) + add_rho_to_equilibrium!(ids) + end + # Check if core_profile is available, if not extrapolate from edge profile + if length(cpp1d) == 0 + fill_in_extrapolated_core_profile(ids, "electrons.density") + elseif length(cpp1d) != length(epggd) + error( + "Number of edge profiles time slices does not match number of " \ + "core profile time slices", + ) + end + fix_eq_time_idx = length(ids.equilibrium.time_slice) == 1 + for ii ∈ eachindex(epggd) + ch.n_e_line.time[ii] = epggd[ii].time + eq_time_idx = fix_eq_time_idx ? 1 : ii + # WIP Integrate electron density along line of sight and add to value + ch.n_e_line.value[ii] + end + end + end +end From 887f882d9e67fe3fcee245defc42f3d27a29b9aa Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Tue, 26 Sep 2023 19:26:25 -0700 Subject: [PATCH 04/12] Synthetic Interferometer calculation added Using QuadGK, a simple line integration of electron density is added. However, no test has been added and this has not been tested yet but the workflow would more or less remain teh same. Also note that for line average density, currently the total path length is used while we want to average by path length in plasma. --- Project.toml | 8 ++++++ src/interferometer.jl | 58 ++++++++++++++++++++++++++++++++++++++----- 2 files changed, 60 insertions(+), 6 deletions(-) diff --git a/Project.toml b/Project.toml index 35acc32..e7d065f 100644 --- a/Project.toml +++ b/Project.toml @@ -2,3 +2,11 @@ name = "SynthDiag" uuid = "c6e34158-aa22-49e4-bb12-15cbcb518ce6" authors = ["Anchal Gupta "] version = "0.1.0" + +[deps] +GGDUtils = "b7b5e640-9b39-4803-84eb-376048795def" +OMAS = "91cfaa06-6526-4804-8666-b540b3feef2f" +PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab" +QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +SD4SOLPS = "d8db6f1b-e564-4c04-bba3-ef399f78c070" +SOLPS2IMAS = "09becab6-0636-4c23-a92a-2b3723265c31" diff --git a/src/interferometer.jl b/src/interferometer.jl index bb00dc2..4363dcc 100644 --- a/src/interferometer.jl +++ b/src/interferometer.jl @@ -1,5 +1,7 @@ import PhysicalConstants.CODATA2018: c_0, ε_0, m_e, e -import SD4SOLPS: fill_in_extrapolated_core_profile, check_rho_1d, add_rho_to_equilibrium +import SD4SOLPS: fill_in_extrapolated_core_profile, check_rho_1d, + add_rho_to_equilibrium, core_profile_2d +import QuadGK: quadgk default_ifo = "$(@__DIR__)/default_interferometer.yml" @@ -37,7 +39,7 @@ Add interferometer to IMAS structure using a Dict and compute the line integrate electron density if not present """ function add_interferometer( - @nospecialize(ids::OMAS.dd)=OMAS.dd(), + @nospecialize(ids::OMAS.dd), config::Dict{Symbol, Any}, )::OMAS.dd OMAS.dict2imas(config, ids) @@ -51,7 +53,7 @@ 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::OMAS.dd)) +function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-6) for ch ∈ ids.interferometer.channel k = [2 * π / ch.wavelength[ii].value for ii ∈ 1:2] for i1 ∈ 1:2 @@ -82,12 +84,56 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd)) "core profile time slices", ) end + # Parametrize line of sight + if ch.line_of_sight.third_point == + OMAS.interferometer__channel___line_of_sight__third_point() + function los(s::Float64) + fp = ch.line_of_sight.first_point.r, ch.line_of_sight.first_point.z + sp = + ch.line_of_sight.second_point.r, ch.line_of_sight.second_point.z + return fp .+ s .* (sp .- fp) + end + else + function los(s::Float64) + fp = ch.line_of_sight.first_point.r, ch.line_of_sight.first_point.z + sp = + ch.line_of_sight.second_point.r, ch.line_of_sight.second_point.z + tp = ch.line_of_sight.third_point.r, ch.line_of_sight.third_point.z + if s <= 0.5 + return fp .+ 2 * s .* (sp .- fp) + else + return sp .+ (s - 0.5) .* (tp .- sp) + end + end + end + # Derivative of length along line of sight with respect to s + dl_ds(s) = + if s <= 0.5 + return sqrt(sum((2 .* (sp .- fp)) .^ 2)) + else + return sqrt(sum((2 .* (sp .- tp)) .^ 2)) + end fix_eq_time_idx = length(ids.equilibrium.time_slice) == 1 for ii ∈ eachindex(epggd) - ch.n_e_line.time[ii] = epggd[ii].time + ch.n_e_line.time[ii] = ch.n_e_line_average.time[ii] = epggd[ii].time eq_time_idx = fix_eq_time_idx ? 1 : ii - # WIP Integrate electron density along line of sight and add to value - ch.n_e_line.value[ii] + integrand(s) = + core_profile_2d( + ids, + ii, + eq_time_idx, + "electrons.density", + los(s)..., + ) * dl_ds.(s) + ch.n_e_line.data[ii] = quadgk(integrand, 0, 1, rtol)[1] + # TO DO: Following line should divide by length in plasma (core) region + ch.n_e_line_average.data[ii] = + ch.n_e_line.data[ii] / quadgk(dl_ds, 0, 1, rtol)[1] + for lam ∈ ch.wavelength + lam.phase_corrected.time[ii] = epggd[ii].time + lam.phase_corrected.data[ii] = + ch.n_e_line.data[ii] / lam.phase_to_n_e_line + end end end end From dc14a0fd12be3e5b6895b4eed39151706575be10 Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Wed, 27 Sep 2023 11:44:00 -0700 Subject: [PATCH 05/12] Adding sample files using dvc --- .dvc/.gitignore | 3 +++ .dvc/config | 0 .dvcignore | 3 +++ sample/.gitignore | 4 ++++ sample/b2fgmtry.dvc | 12 ++++++++++++ sample/b2mn.dat.dvc | 12 ++++++++++++ sample/b2time.nc.dvc | 13 +++++++++++++ sample/g002296.00200.dvc | 12 ++++++++++++ 8 files changed, 59 insertions(+) create mode 100644 .dvc/.gitignore create mode 100644 .dvc/config create mode 100644 .dvcignore create mode 100644 sample/.gitignore create mode 100644 sample/b2fgmtry.dvc create mode 100644 sample/b2mn.dat.dvc create mode 100644 sample/b2time.nc.dvc create mode 100644 sample/g002296.00200.dvc diff --git a/.dvc/.gitignore b/.dvc/.gitignore new file mode 100644 index 0000000..528f30c --- /dev/null +++ b/.dvc/.gitignore @@ -0,0 +1,3 @@ +/config.local +/tmp +/cache diff --git a/.dvc/config b/.dvc/config new file mode 100644 index 0000000..e69de29 diff --git a/.dvcignore b/.dvcignore new file mode 100644 index 0000000..5197305 --- /dev/null +++ b/.dvcignore @@ -0,0 +1,3 @@ +# Add patterns of files dvc should ignore, which could improve +# the performance. Learn more at +# https://dvc.org/doc/user-guide/dvcignore diff --git a/sample/.gitignore b/sample/.gitignore new file mode 100644 index 0000000..0746529 --- /dev/null +++ b/sample/.gitignore @@ -0,0 +1,4 @@ +/b2fgmtry +/b2time.nc +/b2mn.dat +/g002296.00200 diff --git a/sample/b2fgmtry.dvc b/sample/b2fgmtry.dvc new file mode 100644 index 0000000..0701492 --- /dev/null +++ b/sample/b2fgmtry.dvc @@ -0,0 +1,12 @@ +md5: bcf05d7d7b6d91ef1025c0dc018c4bbf +frozen: true +deps: +- path: ITER_Lore_2296_00000/baserun/b2fgmtry + repo: + url: git@github.com:ProjectTorreyPines/SOLPSTestSamples.git + rev_lock: df499f1275428ec06175c48dc4af6ecbd5ec6117 +outs: +- md5: 6589953daeef49b21744d159af970cbb + size: 2908856 + hash: md5 + path: b2fgmtry diff --git a/sample/b2mn.dat.dvc b/sample/b2mn.dat.dvc new file mode 100644 index 0000000..e13ac93 --- /dev/null +++ b/sample/b2mn.dat.dvc @@ -0,0 +1,12 @@ +md5: c0c920771883f03bc0abe98425b38d5e +frozen: true +deps: +- path: ITER_Lore_2296_00000/run_time_dep_EIRENE_jdl_to_ss_cont_sine2_2d_output/b2mn.dat + repo: + url: git@github.com:ProjectTorreyPines/SOLPSTestSamples.git + rev_lock: df499f1275428ec06175c48dc4af6ecbd5ec6117 +outs: +- md5: 1af9d167996ea1a420f49e0d89c98e6d + size: 13925 + hash: md5 + path: b2mn.dat diff --git a/sample/b2time.nc.dvc b/sample/b2time.nc.dvc new file mode 100644 index 0000000..b4a7624 --- /dev/null +++ b/sample/b2time.nc.dvc @@ -0,0 +1,13 @@ +md5: b94876ee58c5d62734d4ee32c28bb0da +frozen: true +deps: +- path: + ITER_Lore_2296_00000/run_time_dep_EIRENE_jdl_to_ss_cont_sine2_2d_output/b2time.nc + repo: + url: git@github.com:ProjectTorreyPines/SOLPSTestSamples.git + rev_lock: df499f1275428ec06175c48dc4af6ecbd5ec6117 +outs: +- md5: 778b3bb62b82f493c785373d872f6a1f + size: 13630424 + hash: md5 + path: b2time.nc diff --git a/sample/g002296.00200.dvc b/sample/g002296.00200.dvc new file mode 100644 index 0000000..7165114 --- /dev/null +++ b/sample/g002296.00200.dvc @@ -0,0 +1,12 @@ +md5: 1d967bf4b5e9d024a87a249859dd10a1 +frozen: true +deps: +- path: ITER_Lore_2296_00000/EQDSK/g002296.00200 + repo: + url: git@github.com:ProjectTorreyPines/SOLPSTestSamples.git + rev_lock: df499f1275428ec06175c48dc4af6ecbd5ec6117 +outs: +- md5: 5ab3f75addd8f1191a49e5243c96d23d + size: 2167720 + hash: md5 + path: g002296.00200 From c2b0e15a45a33a86733e6db522ce7dda08258edf Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Wed, 27 Sep 2023 18:15:23 -0700 Subject: [PATCH 06/12] WIP: Added test, waiting for bug fix in SD4SOLPS * Added a test but it fails due to a bug in SD4SOLPS * Converted default_interferometer.yml to json * Added gridspacedesc.yml in samples * Fixed Project.toml --- Project.toml | 1 + sample/gridspacedesc.yml | 196 ++++++++++++++++++++++++++++++++ src/default_interferometer.json | 122 ++++++++++++++++++++ src/default_interferometer.yml | 82 ------------- src/interferometer.jl | 89 ++++++++------- test/Project.toml | 2 + test/runtests.jl | 13 +++ 7 files changed, 385 insertions(+), 120 deletions(-) create mode 100644 sample/gridspacedesc.yml create mode 100644 src/default_interferometer.json delete mode 100644 src/default_interferometer.yml create mode 100644 test/Project.toml create mode 100644 test/runtests.jl diff --git a/Project.toml b/Project.toml index e7d065f..8a04c9c 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Anchal Gupta "] version = "0.1.0" [deps] +EFIT = "cda752c5-6b03-55a3-9e33-132a441b0c17" GGDUtils = "b7b5e640-9b39-4803-84eb-376048795def" OMAS = "91cfaa06-6526-4804-8666-b540b3feef2f" PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab" diff --git a/sample/gridspacedesc.yml b/sample/gridspacedesc.yml new file mode 100644 index 0000000..0c82f2a --- /dev/null +++ b/sample/gridspacedesc.yml @@ -0,0 +1,196 @@ +identifier: # grid_ggd[#].Identifier properties + name: Sven + index: 1 + description: This is a grid +space: + 1: # space_number + identifier: # grid_ggd[it].space[#].identifier properties + name: linear + index: 1 + description: Linear + geometry_type: # grid_dd[#].space[#].geometry_type properties + name: "standard" + index: 0 # 0: standard, 1:Fourier, >1: Fourier with periodicity + description: "trying to hold a b2/solps mesh here" # Verbose descripition of space + coordinates_type: [4, 3] # r, z +grid_subset: + 1: + dimension: 0 + identifier: + name: nodes + index: 1 + description: All nodes (0D) belonging to the associated spaces, implicit declaration (no need to replicate the grid elements in the grid_subset structure). + 2: + dimension: 1 + identifier: + name: faces + index: 2 + description: All faces (1D) belonging to the associated spaces, implicit declaration (no need to replicate the grid elements in the grid_subset structure) + 3: + dimension: 1 + identifier: + name: x_aligned_faces + index: 3 + description: All x-aligned (polidally aligned) faces + 4: + dimension: 1 + identifier: + name: y_aligned_faces + index: 4 + description: All y-aligned (radially aligned) faces + 5: + dimension: 2 + identifier: + name: cells + index: 5 + description: All cells (2D) belonging to the associated spaces, implicit declaration (no need to replicate the grid elements in the grid_subset structure) + 6: + dimension: 0 + identifier: + name: x_points + index: 6 + description: Nodes defining magnetic X-points (poloidal field nulls) + 7: + dimension: 1 + identifier: + name: core_cut + index: 7 + description: Y-aligned faces (edges) inside the separatrix connecting to the active X-point + 8: + dimension: 1 + identifier: + name: PFR_cut + index: 8 + description: Y-aligned faces (edges) inside the private flux region connecting to the active X-point + 9: + dimension: 1 + identifier: + name: outer_throat + index: 9 + description: Y-aligned faces in the outer SOL connecting to the active X-point + 10: + dimension: 1 + identifier: + name: inner_throat + index: 10 + description: Y-aligned faces in the inner SOL connecting to the active X-point + 11: + dimension: 1 + identifier: + name: outer_midplane + index: 11 + description: Y-aligned faces connecting to the node closest to the outer midplane on the separatrix + 12: + dimension: 1 + identifier: + name: inner_midplane + index: 12 + description: Y-aligned faces connecting to the node closest to the inner midplane on the separatrix + 13: + dimension: 1 + identifier: + name: outer_target + index: 13 + description: Y-aligned faces defining the outer divertor target + 14: + dimension: 1 + identifier: + name: inner_target + index: 14 + description: Y-aligned faces defining the inner divertor target + 15: + dimension: 1 + identifier: + name: core_boundary + index: 15 + description: Innermost X-aligned faces + 16: + dimension: 1 + identifier: + name: separatrix + index: 16 + description: X-aligned faces defining the active separatrix + 17: + dimension: 1 + identifier: + name: main_chamber_wall + index: 17 + description: X-aligned faces defining the main chamber wall outside of the divertor regions + 18: + dimension: 1 + identifier: + name: outer_baffle + index: 18 + description: X-aligned faces defining the chamber wall of the outer active divertor region + 19: + dimension: 1 + identifier: + name: inner_baffle + index: 19 + description: X-aligned faces defining the chamber wall of the inner active divertor region + 20: + dimension: 1 + identifier: + name: outer_PFR_wall + index: 20 + description: X-aligned faces defining the private flux region wall of the outer active divertor region + 21: + dimension: 1 + identifier: + name: inner_PFR_wall + index: 21 + description: X-aligned faces defining the private flux region wall of the inner active divertor region + 22: + dimension: 2 + identifier: + name: core + index: 22 + description: Cells (2D) inside the active separatrix in the same lobe as the magnetic axis + 23: + dimension: 2 + identifier: + name: sol + index: 23 + description: Cells (2D) outside the active separatrix and outside the divertor regions + 24: + dimension: 2 + identifier: + name: outer_divertor + index: 24 + description: Cells (2D) defining the (primary, if there are 2 X-points) outer divertor region + 25: + dimension: 2 + identifier: + name: inner_divertor + index: 25 + description: Cells (2D) defining the (primary, if there are 2 X-points) inner divertor region + 26: + dimension: 1 + identifier: + name: core_sol + index: 26 + description: X-aligned faces defining the part of the active separatrix separating core and SOL + 27: + dimension: 0 + identifier: + name: outer_mid_plane_separatrix + index: 101 + description: Point on the active separatrix at the outer mid-plane + 28: + dimension: 0 + identifier: + name: inner_mid_plane_separatrix + index: 102 + description: Point on the active separatrix at the inner mid-plane + 29: + dimension: 0 + identifier: + name: outer_target_separatrix + index: 103 + description: Point on the active separatrix at the outer active divertor target plate + 30: + dimension: 0 + identifier: + name: inner_target_separatrix + index: 104 + description: Point on the active separatrix at the inner active divertor target plate diff --git a/src/default_interferometer.json b/src/default_interferometer.json new file mode 100644 index 0000000..3f133be --- /dev/null +++ b/src/default_interferometer.json @@ -0,0 +1,122 @@ +{ + "interferometer": { + "channel": [ + { + "name": "V1", + "identifier": "V1", + "line_of_sight": { + "first_point": { + "phi": 0.0, + "r": 5.0, + "z": -5.0 + }, + "second_point": { + "phi": 0.0, + "r": 5.0, + "z": 5.0 + }, + "third_point": { + "phi": 0.0, + "r": 5.0, + "z": -5.0 + } + }, + "wavelength": [ + { + "value": 10.6e-6 + }, + { + "value": 6.33e-7 + } + ] + }, + { + "name": "V2", + "identifier": "V2", + "line_of_sight": { + "first_point": { + "phi": 0.0, + "r": 6.0, + "z": -5.0 + }, + "second_point": { + "phi": 0.0, + "r": 6.0, + "z": 5.0 + }, + "third_point": { + "phi": 0.0, + "r": 6.0, + "z": -5.0 + } + }, + "wavelength": [ + { + "value": 10.6e-6 + }, + { + "value": 6.33e-7 + } + ] + }, + { + "name": "V3", + "identifier": "V3", + "line_of_sight": { + "first_point": { + "phi": 0.0, + "r": 7.0, + "z": -5.0 + }, + "second_point": { + "phi": 0.0, + "r": 7.0, + "z": 5.0 + }, + "third_point": { + "phi": 0.0, + "r": 7.0, + "z": -5.0 + } + }, + "wavelength": [ + { + "value": 10.6e-6 + }, + { + "value": 6.33e-7 + } + ] + }, + { + "name": "H1", + "identifier": "H1", + "line_of_sight": { + "first_point": { + "phi": 0.0, + "r": 9.0, + "z": 0.0 + }, + "second_point": { + "phi": 0.0, + "r": 3.0, + "z": 0.0 + }, + "third_point": { + "phi": 0.0, + "r": 9.0, + "z": 0.0 + } + }, + "wavelength": [ + { + "value": 10.6e-6 + }, + { + "value": 6.33e-7 + } + ] + } + ] + } +} \ No newline at end of file diff --git a/src/default_interferometer.yml b/src/default_interferometer.yml deleted file mode 100644 index 55c0d50..0000000 --- a/src/default_interferometer.yml +++ /dev/null @@ -1,82 +0,0 @@ -interferometer: - channel: - - name: V1 - identifier: 0 - line_of_sight: - first_point: - phi: 0.0 - r: 5.0 - z: -5.0 - second_point: - phi: 0.0 - r: 5.0 - z: 5.0 - third_point: - phi: 0.0 - r: 5.0 - z: -5.0 - wavelength: - - phase_to_n_e_line: - value: 10.6e-6 # in m, for CO2 - - phase_to_n_e_line: - value: 0.633e-6 # in m, for HeNe - - name: V2 - identifier: 1 - line_of_sight: - first_point: - phi: 0.0 - r: 6.0 - z: -5.0 - second_point: - phi: 0.0 - r: 6.0 - z: 5.0 - third_point: - phi: 0.0 - r: 6.0 - z: -5.0 - wavelength: - - phase_to_n_e_line: - value: 10.6e-6 # in m, for CO2 - - phase_to_n_e_line: - value: 0.633e-6 # in m, for HeNe - - name: V3 - identifier: 1 - line_of_sight: - first_point: - phi: 0.0 - r: 7.0 - z: -5.0 - second_point: - phi: 0.0 - r: 7.0 - z: 5.0 - third_point: - phi: 0.0 - r: 7.0 - z: -5.0 - wavelength: - - phase_to_n_e_line: - value: 10.6e-6 # in m, for CO2 - - phase_to_n_e_line: - value: 0.633e-6 # in m, for HeNe - - name: H1 - identifier: 1 - line_of_sight: - first_point: - phi: 0.0 - r: 9.0 - z: 0.0 - second_point: - phi: 0.0 - r: 3.0 - z: 0.0 - third_point: - phi: 0.0 - r: 9.0 - z: 0.0 - wavelength: - - phase_to_n_e_line: - value: 10.6e-6 # in m, for CO2 - - phase_to_n_e_line: - value: 0.633e-6 # in m, for HeNe diff --git a/src/interferometer.jl b/src/interferometer.jl index 4363dcc..58e03bb 100644 --- a/src/interferometer.jl +++ b/src/interferometer.jl @@ -1,9 +1,9 @@ import PhysicalConstants.CODATA2018: c_0, ε_0, m_e, e -import SD4SOLPS: fill_in_extrapolated_core_profile, check_rho_1d, - add_rho_to_equilibrium, core_profile_2d +import SD4SOLPS: fill_in_extrapolated_core_profile!, check_rho_1d, + add_rho_to_equilibrium!, core_profile_2d import QuadGK: quadgk -default_ifo = "$(@__DIR__)/default_interferometer.yml" +default_ifo = "$(@__DIR__)/default_interferometer.json" """ add_interferometer( @@ -12,17 +12,16 @@ default_ifo = "$(@__DIR__)/default_interferometer.yml" )::OMAS.dd -Add interferometer to IMAS structure using either a YAML or JSON file and compute the +Add interferometer to IMAS structure using a JSON file and compute the line integrated electron density if not present """ -function add_interferometer( +function add_interferometer!(config::String=default_ifo, @nospecialize(ids::OMAS.dd)=OMAS.dd(), - config::String=default_ifo, )::OMAS.dd - if endswith(config, ".yml") || endswith(config, ".yaml") - OMAS.yaml2imas(config, ids) - elseif endswith(config, ".json") + if endswith(config, ".json") OMAS.json2imas(config, ids) + else + error("Only JSON files are supported.") end compute_interferometer(ids) return ids @@ -38,9 +37,8 @@ end Add interferometer to IMAS structure using a Dict and compute the line integrated electron density if not present """ -function add_interferometer( - @nospecialize(ids::OMAS.dd), - config::Dict{Symbol, Any}, +function add_interferometer!(config::Dict{Symbol, Any}, + @nospecialize(ids::OMAS.dd)=OMAS.dd(), )::OMAS.dd OMAS.dict2imas(config, ids) compute_interferometer(ids) @@ -62,15 +60,21 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-6) if lam.phase_to_n_e_line == 0 # 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) + ( + 2 * m_e * ε_0 * c_0^2 / e^2 * (k[i1] * k[i2]^2) / + (k[i2]^2 - k[i1]^2) + ).val end end # Special case when measurement has not been made but edge profile data exists if length(ch.n_e_line.time) == 0 && length(ids.edge_profiles.ggd) > 0 epggd = ids.edge_profiles.ggd cpp1d = ids.core_profiles.profiles_1d - ch.n_e_line.time = ch.n_e_line.value = zeros(length(epggd)) + nt = length(epggd) + 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) # Check if equilibrium data contains rho, if not add it if !check_rho_1d(ids; time_slice=1) add_rho_to_equilibrium!(ids) @@ -85,37 +89,31 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-6) ) end # Parametrize line of sight + fp = rzphi2xyz(ch.line_of_sight.first_point) + sp = rzphi2xyz(ch.line_of_sight.second_point) + tp = rzphi2xyz(ch.line_of_sight.third_point) if ch.line_of_sight.third_point == OMAS.interferometer__channel___line_of_sight__third_point() - function los(s::Float64) - fp = ch.line_of_sight.first_point.r, ch.line_of_sight.first_point.z - sp = - ch.line_of_sight.second_point.r, ch.line_of_sight.second_point.z - return fp .+ s .* (sp .- fp) - end + los = (s::Float64) -> xyz2rz((fp .+ s .* (sp .- fp))...) + dl_ds = (s::Float64) -> sqrt(sum((sp .- fp) .^ 2)) else - function los(s::Float64) - fp = ch.line_of_sight.first_point.r, ch.line_of_sight.first_point.z - sp = - ch.line_of_sight.second_point.r, ch.line_of_sight.second_point.z - tp = ch.line_of_sight.third_point.r, ch.line_of_sight.third_point.z - if s <= 0.5 - return fp .+ 2 * s .* (sp .- fp) - else - return sp .+ (s - 0.5) .* (tp .- sp) - end - end - end - # Derivative of length along line of sight with respect to s - dl_ds(s) = - if s <= 0.5 + los = + (s::Float64) -> + if s <= 0.5 + return xyz2rz((fp .+ 2 * s .* (sp .- fp))...) + else + return xyz2rz((sp .+ (s - 0.5) .* (tp .- sp))...) + end + dl_ds = (s::Float64) -> if s <= 0.5 return sqrt(sum((2 .* (sp .- fp)) .^ 2)) else return sqrt(sum((2 .* (sp .- tp)) .^ 2)) end + end fix_eq_time_idx = length(ids.equilibrium.time_slice) == 1 for ii ∈ eachindex(epggd) - ch.n_e_line.time[ii] = ch.n_e_line_average.time[ii] = epggd[ii].time + ch.n_e_line.time[ii] = epggd[ii].time + ch.n_e_line_average.time[ii] = epggd[ii].time eq_time_idx = fix_eq_time_idx ? 1 : ii integrand(s) = core_profile_2d( @@ -124,7 +122,8 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-6) eq_time_idx, "electrons.density", los(s)..., - ) * dl_ds.(s) + ) .* dl_ds(s) + println(integrand(0.5)) ch.n_e_line.data[ii] = quadgk(integrand, 0, 1, rtol)[1] # TO DO: Following line should divide by length in plasma (core) region ch.n_e_line_average.data[ii] = @@ -138,3 +137,17 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-6) end end end + +function rzphi2xyz( + point::Union{OMAS.interferometer__channel___line_of_sight__first_point, + OMAS.interferometer__channel___line_of_sight__second_point, + OMAS.interferometer__channel___line_of_sight__third_point}, +) + r, z, phi = point.r, point.z, point.phi + return r * cos(phi), r * sin(phi), z +end + +function xyz2rz(x::Float64, y::Float64, z::Float64) + r = sqrt(x^2 + y^2) + return r, z +end diff --git a/test/Project.toml b/test/Project.toml new file mode 100644 index 0000000..0c36332 --- /dev/null +++ b/test/Project.toml @@ -0,0 +1,2 @@ +[deps] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..0baf7af --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,13 @@ +using SynthDiag: add_interferometer! +using SD4SOLPS: preparation +using Test + +@testset "interferometer" begin + eqdsk_file = "g002296.00200" + sample_dir = "$(@__DIR__)/../sample" + test_dir = mktempdir() + ids = preparation(eqdsk_file, sample_dir; filename=test_dir * "/output") + add_interferometer!("$(@__DIR__)/../src/default_interferometer.json", ids) + # Just checking if the function runs through for now + @test true +end From 7a1e234a0ec857428de2685b7d6a83846ebc1f40 Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Fri, 29 Sep 2023 18:21:42 -0700 Subject: [PATCH 07/12] Using diff extrapolation based on loc, works, slow MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Utilizing the new ∈ operator from GGDUtils.jl, now using core_profile_2d for getting values inside core region and using GGDUtils.interp for getting values inside SOLPS grid outside of core region. Everywhere else, teh electron density is assumed to be zero. Chord length inside core is also measured using the same functionality. However, the implementation is very slow. It takes about 4s to integrate electron density along a chord for each time step . This is not acceptable. I'll look deeper into what is causing the slow down and fix it. --- src/interferometer.jl | 92 +++++++++++++++++++++++++++++++++++++++---- 1 file changed, 85 insertions(+), 7 deletions(-) diff --git a/src/interferometer.jl b/src/interferometer.jl index 58e03bb..a07c341 100644 --- a/src/interferometer.jl +++ b/src/interferometer.jl @@ -2,6 +2,8 @@ import PhysicalConstants.CODATA2018: c_0, ε_0, m_e, e import SD4SOLPS: fill_in_extrapolated_core_profile!, check_rho_1d, add_rho_to_equilibrium!, core_profile_2d import QuadGK: quadgk +using GGDUtils +using SOLPS2IMAS: SOLPS2IMAS default_ifo = "$(@__DIR__)/default_interferometer.json" @@ -51,7 +53,7 @@ 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::OMAS.dd), rtol::Float64=1e-6) +function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1.0) for ch ∈ ids.interferometer.channel k = [2 * π / ch.wavelength[ii].value for ii ∈ 1:2] for i1 ∈ 1:2 @@ -75,6 +77,10 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-6) 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 # Check if equilibrium data contains rho, if not add it if !check_rho_1d(ids; time_slice=1) add_rho_to_equilibrium!(ids) @@ -111,23 +117,55 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-6) end end fix_eq_time_idx = length(ids.equilibrium.time_slice) == 1 + fix_ep_grid_ggd_idx = length(ids.edge_profiles.grid_ggd) == 1 for ii ∈ eachindex(epggd) ch.n_e_line.time[ii] = epggd[ii].time ch.n_e_line_average.time[ii] = epggd[ii].time eq_time_idx = fix_eq_time_idx ? 1 : ii + ep_grid_ggd_idx = fix_ep_grid_ggd_idx ? 1 : ii + ep_grid_ggd = ids.edge_profiles.grid_ggd[ep_grid_ggd_idx] + ep_space = ep_grid_ggd.space[1] + cells = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 5) + core = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 22) + sol = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 23) + + sep_core = OMAS.edge_profiles__grid_ggd___grid_subset() + sep_core.element = + SOLPS2IMAS.subset_do( + intersect, + SOLPS2IMAS.get_subset_boundary(ep_space, sol), + SOLPS2IMAS.get_subset_boundary(ep_space, core), + ) + + SOLPS_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() + SOLPS_bnd.element = SOLPS2IMAS.get_subset_boundary(ep_space, cells) + ep_kdtree = get_kdtree(ep_space) + ep_n_e = interp( + get_prop_with_grid_subset_index(epggd[ii].electrons.density, 5), + ep_kdtree, + ) integrand(s) = - core_profile_2d( + n_e( ids, ii, eq_time_idx, - "electrons.density", + sep_core, + SOLPS_bnd, + ep_space, + ep_n_e, los(s)..., ) .* dl_ds(s) - println(integrand(0.5)) - ch.n_e_line.data[ii] = quadgk(integrand, 0, 1, rtol)[1] - # TO DO: Following line should divide by length in plasma (core) region + println("Calculating integrated electron density") + @time ch.n_e_line.data[ii] = quadgk(integrand, 0, 1, rtol)[1] + println(ch.n_e_line.data[ii]) + core_chord_integrand(s) = + is_in_core(sep_core, ep_space, los(s)...) .* dl_ds(s) + println("Calculating core chord length") + @time core_chord_length = quadgk(core_chord_integrand, 0, 1, rtol)[1] + println(core_chord_length) ch.n_e_line_average.data[ii] = - ch.n_e_line.data[ii] / quadgk(dl_ds, 0, 1, rtol)[1] + ch.n_e_line.data[ii] / core_chord_length + println("Average: ", ch.n_e_line_average.data[ii]) for lam ∈ ch.wavelength lam.phase_corrected.time[ii] = epggd[ii].time lam.phase_corrected.data[ii] = @@ -151,3 +189,43 @@ function xyz2rz(x::Float64, y::Float64, z::Float64) r = sqrt(x^2 + y^2) return r, z end + +function n_e( + @nospecialize(ids::OMAS.dd), + prof_time_idx::Int64, + eq_time_idx::Int64, + sep_core::OMAS.edge_profiles__grid_ggd___grid_subset, + SOLPS_bnd::OMAS.edge_profiles__grid_ggd___grid_subset, + ep_space::OMAS.edge_profiles__grid_ggd___space, + ep_n_e, + r::Float64, + z::Float64, +) + if (r, z) ∈ (sep_core, ep_space) + return core_profile_2d( + ids, + prof_time_idx, + eq_time_idx, + "electrons.density", + r, + z, + ) + elseif (r, z) ∈ (SOLPS_bnd, ep_space) + return ep_n_e(r, z) + else + return 0.0 + end +end + +function is_in_core( + sep_core::OMAS.edge_profiles__grid_ggd___grid_subset, + ep_space::OMAS.edge_profiles__grid_ggd___space, + r::Float64, + z::Float64, +) + if (r, z) ∈ (sep_core, ep_space) + return 1.0 + else + return 0.0 + end +end \ No newline at end of file From 5415e58dbf17423bd46175d3d07b3815cec44e7d Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Mon, 2 Oct 2023 17:46:38 -0700 Subject: [PATCH 08/12] Using exported core_profile_2d function --- src/interferometer.jl | 37 ++++++++++++++----------------------- 1 file changed, 14 insertions(+), 23 deletions(-) diff --git a/src/interferometer.jl b/src/interferometer.jl index a07c341..48bdc22 100644 --- a/src/interferometer.jl +++ b/src/interferometer.jl @@ -53,7 +53,7 @@ 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::OMAS.dd), rtol::Float64=1.0) +function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-6) for ch ∈ ids.interferometer.channel k = [2 * π / ch.wavelength[ii].value for ii ∈ 1:2] for i1 ∈ 1:2 @@ -125,7 +125,10 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1.0) ep_grid_ggd_idx = fix_ep_grid_ggd_idx ? 1 : ii ep_grid_ggd = ids.edge_profiles.grid_ggd[ep_grid_ggd_idx] ep_space = ep_grid_ggd.space[1] - cells = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 5) + all_cells = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 5) + SOLPS_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() + SOLPS_bnd.element = SOLPS2IMAS.get_subset_boundary(ep_space, all_cells) + core_bnd = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 15) core = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 22) sol = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 23) @@ -137,22 +140,19 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1.0) SOLPS2IMAS.get_subset_boundary(ep_space, core), ) - SOLPS_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() - SOLPS_bnd.element = SOLPS2IMAS.get_subset_boundary(ep_space, cells) ep_kdtree = get_kdtree(ep_space) ep_n_e = interp( get_prop_with_grid_subset_index(epggd[ii].electrons.density, 5), ep_kdtree, ) + cp_n_e = core_profile_2d(ids, ii, eq_time_idx, "electrons.density") integrand(s) = n_e( - ids, - ii, - eq_time_idx, - sep_core, + cp_n_e, + ep_n_e, + core_bnd, SOLPS_bnd, ep_space, - ep_n_e, los(s)..., ) .* dl_ds(s) println("Calculating integrated electron density") @@ -191,25 +191,16 @@ function xyz2rz(x::Float64, y::Float64, z::Float64) end function n_e( - @nospecialize(ids::OMAS.dd), - prof_time_idx::Int64, - eq_time_idx::Int64, - sep_core::OMAS.edge_profiles__grid_ggd___grid_subset, + cp_n_e::Function, + ep_n_e::Function, + core_bnd::OMAS.edge_profiles__grid_ggd___grid_subset, SOLPS_bnd::OMAS.edge_profiles__grid_ggd___grid_subset, ep_space::OMAS.edge_profiles__grid_ggd___space, - ep_n_e, r::Float64, z::Float64, ) - if (r, z) ∈ (sep_core, ep_space) - return core_profile_2d( - ids, - prof_time_idx, - eq_time_idx, - "electrons.density", - r, - z, - ) + if (r, z) ∈ (core_bnd, ep_space) + return cp_n_e(r, z) elseif (r, z) ∈ (SOLPS_bnd, ep_space) return ep_n_e(r, z) else From 0d6a4443b9ac36581dcee26655f7da85b520668b Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Wed, 4 Oct 2023 11:30:23 -0700 Subject: [PATCH 09/12] Using prepacked core profile interpolation func Using the new core profile interpolation function that returns an inteprolating function instead of values at specific points. This helped increase the speed of the code such that for relative tolerance of 10 minutes, most integrations finish withing 10 seconds. However, the horizontal line of sight chord is taking longer. The major reason for slow speed of numerical integration is the presence of peaks that require a lot of fine resolution function calls by quadGK to reach to the desired relative tolerance. These peaks are naturally present in the data and we can not do much about it. There is bit of commented code in this commit where I tried to use multi processing using threads but that is returning negative inegrated electron density and is not really speeding up the performance, so I have commented it out. Maybe in future this can be done more efficiently --- src/interferometer.jl | 126 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 109 insertions(+), 17 deletions(-) diff --git a/src/interferometer.jl b/src/interferometer.jl index 48bdc22..0b78893 100644 --- a/src/interferometer.jl +++ b/src/interferometer.jl @@ -1,7 +1,7 @@ import PhysicalConstants.CODATA2018: c_0, ε_0, m_e, e import SD4SOLPS: fill_in_extrapolated_core_profile!, check_rho_1d, add_rho_to_equilibrium!, core_profile_2d -import QuadGK: quadgk +import QuadGK: quadgk, BatchIntegrand using GGDUtils using SOLPS2IMAS: SOLPS2IMAS @@ -53,7 +53,7 @@ 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::OMAS.dd), rtol::Float64=1e-6) +function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-1) for ch ∈ ids.interferometer.channel k = [2 * π / ch.wavelength[ii].value for ii ∈ 1:2] for i1 ∈ 1:2 @@ -102,19 +102,21 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-6) OMAS.interferometer__channel___line_of_sight__third_point() los = (s::Float64) -> xyz2rz((fp .+ s .* (sp .- fp))...) dl_ds = (s::Float64) -> sqrt(sum((sp .- fp) .^ 2)) + chord_points = (fp, sp) else los = (s::Float64) -> if s <= 0.5 return xyz2rz((fp .+ 2 * s .* (sp .- fp))...) else - return xyz2rz((sp .+ (s - 0.5) .* (tp .- sp))...) + return xyz2rz((sp .+ 2 * (s - 0.5) .* (tp .- sp))...) end dl_ds = (s::Float64) -> if s <= 0.5 return sqrt(sum((2 .* (sp .- fp)) .^ 2)) else return sqrt(sum((2 .* (sp .- tp)) .^ 2)) end + chord_points = (fp, sp, tp) end fix_eq_time_idx = length(ids.equilibrium.time_slice) == 1 fix_ep_grid_ggd_idx = length(ids.edge_profiles.grid_ggd) == 1 @@ -140,6 +142,11 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-6) SOLPS2IMAS.get_subset_boundary(ep_space, core), ) + valid_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() + valid_bnd.element = + SOLPS2IMAS.subset_do(setdiff, SOLPS_bnd.element, core_bnd.element) + valid_segs = get_intersections(valid_bnd, ep_space, chord_points, los) + println(valid_segs) ep_kdtree = get_kdtree(ep_space) ep_n_e = interp( get_prop_with_grid_subset_index(epggd[ii].electrons.density, 5), @@ -155,13 +162,35 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-6) ep_space, los(s)..., ) .* dl_ds(s) + function integrand!(y, s) + n = Threads.nthreads() + Threads.@threads for i ∈ 1:n + y[i:n:end] .= integrand.(@view(s[i:n:end])) + end + end println("Calculating integrated electron density") + ch.n_e_line.data[ii] = 0.0 + # for seg ∈ valid_segs + # @time ch.n_e_line.data[ii] += + # quadgk( + # BatchIntegrand{Float64}(integrand!), + # seg[1], + # seg[2], + # rtol, + # )[1] + # end @time ch.n_e_line.data[ii] = quadgk(integrand, 0, 1, rtol)[1] println(ch.n_e_line.data[ii]) - core_chord_integrand(s) = - is_in_core(sep_core, ep_space, los(s)...) .* dl_ds(s) - println("Calculating core chord length") - @time core_chord_length = quadgk(core_chord_integrand, 0, 1, rtol)[1] + chord_in_core = get_intersections(sep_core, ep_space, chord_points, los) + core_chord_length = 0.0 + for seg ∈ chord_in_core + core_chord_length += + (seg[2] - seg[1]) * dl_ds((seg[1] + seg[2]) / 2) + end + # core_chord_integrand(s) = + # is_in_core(sep_core, ep_space, los(s)...) .* dl_ds(s) + # println("Calculating core chord length") + # @time core_chord_length = quadgk(core_chord_integrand, 0, 1, rtol)[1] println(core_chord_length) ch.n_e_line_average.data[ii] = ch.n_e_line.data[ii] / core_chord_length @@ -208,15 +237,78 @@ function n_e( end end -function is_in_core( - sep_core::OMAS.edge_profiles__grid_ggd___grid_subset, - ep_space::OMAS.edge_profiles__grid_ggd___space, - r::Float64, - z::Float64, -) - if (r, z) ∈ (sep_core, ep_space) - return 1.0 +function get_intersections(subset, space, points, los) + nodes = space.objects_per_dimension[1].object + edges = space.objects_per_dimension[2].object + if length(points) == 2 + fp, sp = points + three_points = false else - return 0.0 + fp, sp, tp = points + three_points = true + end + segs = [0.0] + for ele ∈ subset.element + edge = edges[ele.object[1].index] + e1 = Tuple(nodes[edge.nodes[1]].geometry) + e2 = Tuple(nodes[edge.nodes[2]].geometry) + s, s2 = intersection_s(e1, e2, fp, sp) + if 0 <= s < 1 && 0 <= s2 < 1 + if three_points + append!(segs, s / 2) + s, s2 = intersection_s(e1, e2, sp, tp) + if 0 <= s < 1 && 0 <= s2 < 1 + append!(segs, s / 2 + 0.5) + end + else + append!(segs, s) + end + end end -end \ No newline at end of file + append!(segs, 1.0) + sort!(segs) + filt_segs = Tuple{Float64, Float64}[] + for ii ∈ 1:length(segs)-1 + check_at = (segs[ii] + segs[ii+1]) / 2 + # print("Seg : ", los(segs[ii]), " to ", los(segs[ii+1]), ", Mid: ") + # print(check_at, ": ", los(check_at), ": ") + if los(check_at) ∈ (subset, space) + append!(filt_segs, [(segs[ii], segs[ii+1])]) + # println("Present") + # else + # println("Absent") + end + end + return filt_segs +end + +function intersection_s( + fp::Tuple{Float64, Float64}, + sp::Tuple{Float64, Float64}, + l1::Tuple{Float64, Float64}, + l2::Tuple{Float64, Float64}, +) + den1 = (sp[1] - fp[1]) * (l2[2] - l1[2]) - (sp[2] - fp[2]) * (l2[1] - l1[1]) + num1 = (sp[1] - fp[1]) * (fp[2] - l1[2]) - (sp[2] - fp[2]) * (fp[1] - l1[1]) + den2 = (l2[1] - l1[1]) * (sp[2] - fp[2]) - (l2[2] - l1[2]) * (sp[1] - fp[1]) + num2 = (l2[1] - l1[1]) * (l1[2] - fp[2]) - (l2[2] - l1[2]) * (l1[1] - fp[1]) + return num1 / den1, num2 / den2 +end + +function intersection_s( + fp::Tuple{Float64, Float64, Float64}, + sp::Tuple{Float64, Float64, Float64}, + l1::Tuple{Float64, Float64}, + l2::Tuple{Float64, Float64}, +) + return intersection_s(xyz2rz(fp...), xyz2rz(sp...), l1, l2) +end + +function intersection_s( + fp::Tuple{Float64, Float64}, + sp::Tuple{Float64, Float64}, + l1::Tuple{Float64, Float64, Float64}, + l2::Tuple{Float64, Float64, Float64}, +) + return intersection_s(fp, sp, xyz2rz(l1...), xyz2rz(l2...)) +end From b339a7d29e618dd1e8b65aba6fd59a0ef1fb3fce Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Mon, 9 Oct 2023 17:07:29 -0700 Subject: [PATCH 10/12] Interferometer now functioning fast enough MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Most changes were required in GGDUtils for making interpolations and projection faster. * The main change in interferometer.jl was to reuse the calculations as much as possible. * I also tried calculating the integrated electron density in each segment of chord (where each segment is either in edge profile or in core profile) but this was slower than the current method of simply integrating over the whole chord and using ∈ to find the relevant function to use for interpolation. * I have left the other method as commented out code right now. Will clean up in next commit to store this work as well. --- src/interferometer.jl | 175 +++++++++++++++++++++++++++++++----------- 1 file changed, 130 insertions(+), 45 deletions(-) diff --git a/src/interferometer.jl b/src/interferometer.jl index 0b78893..0a9fa98 100644 --- a/src/interferometer.jl +++ b/src/interferometer.jl @@ -4,6 +4,7 @@ import SD4SOLPS: fill_in_extrapolated_core_profile!, check_rho_1d, import QuadGK: quadgk, BatchIntegrand using GGDUtils using SOLPS2IMAS: SOLPS2IMAS +using ProgressBars: ProgressBar default_ifo = "$(@__DIR__)/default_interferometer.json" @@ -54,6 +55,37 @@ end - Compute the line integrated electron density if not present """ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-1) + 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] + all_cells = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 5) + SOLPS_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() + SOLPS_bnd.element = SOLPS2IMAS.get_subset_boundary(ep_space, all_cells) + core_bnd = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 15) + core = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 22) + sol = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 23) + + cp_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() + cp_bnd.element = + SOLPS2IMAS.subset_do( + intersect, + SOLPS2IMAS.get_subset_boundary(ep_space, sol), + SOLPS2IMAS.get_subset_boundary(ep_space, core), + ) + SOLPS_out_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() + SOLPS_out_bnd.element = + SOLPS2IMAS.subset_do(setdiff, SOLPS_bnd.element, core_bnd.element) + ep_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() + ep_bnd.element = + SOLPS2IMAS.subset_do(union, SOLPS_out_bnd.element, cp_bnd.element) + + TPS_mats_all_cells = get_TPS_mats(ep_grid_ggd, 5) + + ep_n_e_list = [] + cp_n_e_list = [] + for ch ∈ ids.interferometer.channel k = [2 * π / ch.wavelength[ii].value for ii ∈ 1:2] for i1 ∈ 1:2 @@ -118,41 +150,94 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-1) end chord_points = (fp, sp, tp) end - fix_eq_time_idx = length(ids.equilibrium.time_slice) == 1 - fix_ep_grid_ggd_idx = length(ids.edge_profiles.grid_ggd) == 1 + + # ep_segs = get_intersections(ep_bnd, ep_space, chord_points, los) + # cp_segs = get_intersections(cp_bnd, ep_space, chord_points, los) + for ii ∈ eachindex(epggd) ch.n_e_line.time[ii] = epggd[ii].time ch.n_e_line_average.time[ii] = epggd[ii].time eq_time_idx = fix_eq_time_idx ? 1 : ii - ep_grid_ggd_idx = fix_ep_grid_ggd_idx ? 1 : ii - ep_grid_ggd = ids.edge_profiles.grid_ggd[ep_grid_ggd_idx] - ep_space = ep_grid_ggd.space[1] - all_cells = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 5) - SOLPS_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() - SOLPS_bnd.element = SOLPS2IMAS.get_subset_boundary(ep_space, all_cells) - core_bnd = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 15) - core = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 22) - sol = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 23) + if !fix_ep_grid_ggd_idx + ep_grid_ggd = ids.edge_profiles.grid_ggd[ii] + ep_space = ep_grid_ggd.space[1] + all_cells = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 5) + SOLPS_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() + SOLPS_bnd.element = + SOLPS2IMAS.get_subset_boundary(ep_space, all_cells) + core_bnd = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 15) + core = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 22) + sol = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 23) + + cp_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() + cp_bnd.element = + SOLPS2IMAS.subset_do( + intersect, + SOLPS2IMAS.get_subset_boundary(ep_space, sol), + SOLPS2IMAS.get_subset_boundary(ep_space, core), + ) + SOLPS_out_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() + SOLPS_out_bnd.element = + SOLPS2IMAS.subset_do( + setdiff, + SOLPS_bnd.element, + core_bnd.element, + ) + ep_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() + ep_bnd.element = + SOLPS2IMAS.subset_do( + union, + SOLPS_out_bnd.element, + cp_bnd.element, + ) + + TPS_mats_all_cells = get_TPS_mats(ep_grid_ggd, 5) + + ep_segs = get_intersections(ep_bnd, ep_space, chord_points, los) + cp_segs = get_intersections(cp_bnd, ep_space, chord_points, los) + end - sep_core = OMAS.edge_profiles__grid_ggd___grid_subset() - sep_core.element = - SOLPS2IMAS.subset_do( - intersect, - SOLPS2IMAS.get_subset_boundary(ep_space, sol), - SOLPS2IMAS.get_subset_boundary(ep_space, core), + if length(ep_n_e_list) < ii + append!( + ep_n_e_list, + [interp(epggd[ii].electrons.density, TPS_mats_all_cells, 5)], ) + end + if length(cp_n_e_list) < ii + cp_p1d = ids.core_profiles.profiles_1d[ii] + eqt = ids.equilibrium.time_slice[eq_time_idx] + append!( + cp_n_e_list, + [interp(cp_p1d.electrons.density, cp_p1d, eqt)], + ) + end - valid_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() - valid_bnd.element = - SOLPS2IMAS.subset_do(setdiff, SOLPS_bnd.element, core_bnd.element) - valid_segs = get_intersections(valid_bnd, ep_space, chord_points, los) - println(valid_segs) - ep_kdtree = get_kdtree(ep_space) - ep_n_e = interp( - get_prop_with_grid_subset_index(epggd[ii].electrons.density, 5), - ep_kdtree, - ) - cp_n_e = core_profile_2d(ids, ii, eq_time_idx, "electrons.density") + ep_n_e = ep_n_e_list[ii] + cp_n_e = cp_n_e_list[ii] + + # integ_ep = (s) -> ep_n_e(los(s)...) * dl_ds(s) + # integ_cp = (s) -> cp_n_e.(los(s)...) * dl_ds(s) + + # int_n_e = 0 + # # println("In Edge Profiles:") + # for seg ∈ ep_segs + # # println(seg[1], " to ", seg[2]) + # # println("Slope: ", dl_ds((seg[1] + seg[2]) / 2)) + # # println("Inegrand: ", integ_ep((seg[1] + seg[2]) / 2)) + # seg_int = quadgk(integ_ep, seg[1], seg[2]; rtol=rtol)[1] + # int_n_e += seg_int + # # println("This segment: ", seg_int, "Total: ", int_n_e) + # end + # # println("In Core Profiles:") + # for seg ∈ cp_segs + # # println(seg[1], " to ", seg[2]) + # # println("Slope: ", dl_ds((seg[1] + seg[2]) / 2)) + # # println("Inegrand: ", integ_cp((seg[1] + seg[2]) / 2)) + # seg_int = quadgk(integ_cp, seg[1], seg[2]; rtol=rtol)[1] + # int_n_e += seg_int + # # println("This segment: ", seg_int, "Total: ", int_n_e) + # end + # ch.n_e_line.data[ii] = int_n_e integrand(s) = n_e( cp_n_e, @@ -162,14 +247,14 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-1) ep_space, los(s)..., ) .* dl_ds(s) - function integrand!(y, s) - n = Threads.nthreads() - Threads.@threads for i ∈ 1:n - y[i:n:end] .= integrand.(@view(s[i:n:end])) - end - end - println("Calculating integrated electron density") - ch.n_e_line.data[ii] = 0.0 + # function integrand!(y, s) + # n = Threads.nthreads() + # Threads.@threads for i ∈ 1:n + # y[i:n:end] .= integrand.(@view(s[i:n:end])) + # end + # end + # println("Calculating integrated electron density") + # ch.n_e_line.data[ii] = 0.0 # for seg ∈ valid_segs # @time ch.n_e_line.data[ii] += # quadgk( @@ -179,22 +264,22 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-1) # rtol, # )[1] # end - @time ch.n_e_line.data[ii] = quadgk(integrand, 0, 1, rtol)[1] - println(ch.n_e_line.data[ii]) - chord_in_core = get_intersections(sep_core, ep_space, chord_points, los) + ch.n_e_line.data[ii] = quadgk(integrand, 0, 1; rtol=rtol)[1] + # println(ch.n_e_line.data[ii]) + chord_in_core = get_intersections(cp_bnd, ep_space, chord_points, los) core_chord_length = 0.0 for seg ∈ chord_in_core core_chord_length += (seg[2] - seg[1]) * dl_ds((seg[1] + seg[2]) / 2) end # core_chord_integrand(s) = - # is_in_core(sep_core, ep_space, los(s)...) .* dl_ds(s) + # is_in_core( cp_bnd, ep_space, los(s)...) .* dl_ds(s) # println("Calculating core chord length") # @time core_chord_length = quadgk(core_chord_integrand, 0, 1, rtol)[1] - println(core_chord_length) + # println(core_chord_length) ch.n_e_line_average.data[ii] = ch.n_e_line.data[ii] / core_chord_length - println("Average: ", ch.n_e_line_average.data[ii]) + # println("Average: ", ch.n_e_line_average.data[ii]) for lam ∈ ch.wavelength lam.phase_corrected.time[ii] = epggd[ii].time lam.phase_corrected.data[ii] = @@ -230,10 +315,10 @@ function n_e( ) if (r, z) ∈ (core_bnd, ep_space) return cp_n_e(r, z) - elseif (r, z) ∈ (SOLPS_bnd, ep_space) + else# if (r, z) ∈ (SOLPS_bnd, ep_space) return ep_n_e(r, z) - else - return 0.0 + # else + # return 0.0 end end From 924dbae1fe511f7689e7ff1297833fd39fde3c59 Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Mon, 9 Oct 2023 18:14:24 -0700 Subject: [PATCH 11/12] Cleaned up, same logic and code --- src/interferometer.jl | 231 +++++++++++++----------------------------- test/runtests.jl | 19 ++++ 2 files changed, 91 insertions(+), 159 deletions(-) diff --git a/src/interferometer.jl b/src/interferometer.jl index 0a9fa98..9c33833 100644 --- a/src/interferometer.jl +++ b/src/interferometer.jl @@ -4,7 +4,6 @@ import SD4SOLPS: fill_in_extrapolated_core_profile!, check_rho_1d, import QuadGK: quadgk, BatchIntegrand using GGDUtils using SOLPS2IMAS: SOLPS2IMAS -using ProgressBars: ProgressBar default_ifo = "$(@__DIR__)/default_interferometer.json" @@ -18,7 +17,8 @@ default_ifo = "$(@__DIR__)/default_interferometer.json" Add interferometer to IMAS structure using a JSON file and compute the line integrated electron density if not present """ -function add_interferometer!(config::String=default_ifo, +function add_interferometer!( + config::String=default_ifo, @nospecialize(ids::OMAS.dd)=OMAS.dd(), )::OMAS.dd if endswith(config, ".json") @@ -40,7 +40,8 @@ end Add interferometer to IMAS structure using a Dict and compute the line integrated electron density if not present """ -function add_interferometer!(config::Dict{Symbol, Any}, +function add_interferometer!( + config::Dict{Symbol, Any}, @nospecialize(ids::OMAS.dd)=OMAS.dd(), )::OMAS.dd OMAS.dict2imas(config, ids) @@ -54,33 +55,13 @@ 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::OMAS.dd), rtol::Float64=1e-1) +function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-3) 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] - all_cells = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 5) - SOLPS_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() - SOLPS_bnd.element = SOLPS2IMAS.get_subset_boundary(ep_space, all_cells) - core_bnd = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 15) - core = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 22) - sol = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 23) - - cp_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() - cp_bnd.element = - SOLPS2IMAS.subset_do( - intersect, - SOLPS2IMAS.get_subset_boundary(ep_space, sol), - SOLPS2IMAS.get_subset_boundary(ep_space, core), - ) - SOLPS_out_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() - SOLPS_out_bnd.element = - SOLPS2IMAS.subset_do(setdiff, SOLPS_bnd.element, core_bnd.element) - ep_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() - ep_bnd.element = - SOLPS2IMAS.subset_do(union, SOLPS_out_bnd.element, cp_bnd.element) - + sep_bnd, core_bnd = get_sep_core_bnd(ep_grid_ggd) TPS_mats_all_cells = get_TPS_mats(ep_grid_ggd, 5) ep_n_e_list = [] @@ -132,69 +113,25 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-1) tp = rzphi2xyz(ch.line_of_sight.third_point) if ch.line_of_sight.third_point == OMAS.interferometer__channel___line_of_sight__third_point() - los = (s::Float64) -> xyz2rz((fp .+ s .* (sp .- fp))...) - dl_ds = (s::Float64) -> sqrt(sum((sp .- fp) .^ 2)) chord_points = (fp, sp) else - los = - (s::Float64) -> - if s <= 0.5 - return xyz2rz((fp .+ 2 * s .* (sp .- fp))...) - else - return xyz2rz((sp .+ 2 * (s - 0.5) .* (tp .- sp))...) - end - dl_ds = (s::Float64) -> if s <= 0.5 - return sqrt(sum((2 .* (sp .- fp)) .^ 2)) - else - return sqrt(sum((2 .* (sp .- tp)) .^ 2)) - end chord_points = (fp, sp, tp) end - - # ep_segs = get_intersections(ep_bnd, ep_space, chord_points, los) - # cp_segs = get_intersections(cp_bnd, ep_space, chord_points, los) + los, dl_ds = get_line_of_sight(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 ch.n_e_line_average.time[ii] = epggd[ii].time eq_time_idx = fix_eq_time_idx ? 1 : ii if !fix_ep_grid_ggd_idx + # If grid_ggd is evolving with time, update boundaries ep_grid_ggd = ids.edge_profiles.grid_ggd[ii] ep_space = ep_grid_ggd.space[1] - all_cells = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 5) - SOLPS_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() - SOLPS_bnd.element = - SOLPS2IMAS.get_subset_boundary(ep_space, all_cells) - core_bnd = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 15) - core = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 22) - sol = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 23) - - cp_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() - cp_bnd.element = - SOLPS2IMAS.subset_do( - intersect, - SOLPS2IMAS.get_subset_boundary(ep_space, sol), - SOLPS2IMAS.get_subset_boundary(ep_space, core), - ) - SOLPS_out_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() - SOLPS_out_bnd.element = - SOLPS2IMAS.subset_do( - setdiff, - SOLPS_bnd.element, - core_bnd.element, - ) - ep_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() - ep_bnd.element = - SOLPS2IMAS.subset_do( - union, - SOLPS_out_bnd.element, - cp_bnd.element, - ) - + sep_bnd, core_bnd = get_sep_core_bnd(ep_grid_ggd) TPS_mats_all_cells = get_TPS_mats(ep_grid_ggd, 5) - - ep_segs = get_intersections(ep_bnd, ep_space, chord_points, los) - cp_segs = get_intersections(cp_bnd, ep_space, chord_points, los) + core_chord_length = + get_core_chord_length(sep_bnd, ep_space, chord_points) end if length(ep_n_e_list) < ii @@ -215,71 +152,17 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-1) ep_n_e = ep_n_e_list[ii] cp_n_e = cp_n_e_list[ii] - # integ_ep = (s) -> ep_n_e(los(s)...) * dl_ds(s) - # integ_cp = (s) -> cp_n_e.(los(s)...) * dl_ds(s) + function integrand(s::Real) + r, z = los(s) + if (r, z) ∈ (sep_bnd, ep_space) + return cp_n_e(r, z) .* dl_ds(s) + else + return ep_n_e(r, z) .* dl_ds(s) + end + end - # int_n_e = 0 - # # println("In Edge Profiles:") - # for seg ∈ ep_segs - # # println(seg[1], " to ", seg[2]) - # # println("Slope: ", dl_ds((seg[1] + seg[2]) / 2)) - # # println("Inegrand: ", integ_ep((seg[1] + seg[2]) / 2)) - # seg_int = quadgk(integ_ep, seg[1], seg[2]; rtol=rtol)[1] - # int_n_e += seg_int - # # println("This segment: ", seg_int, "Total: ", int_n_e) - # end - # # println("In Core Profiles:") - # for seg ∈ cp_segs - # # println(seg[1], " to ", seg[2]) - # # println("Slope: ", dl_ds((seg[1] + seg[2]) / 2)) - # # println("Inegrand: ", integ_cp((seg[1] + seg[2]) / 2)) - # seg_int = quadgk(integ_cp, seg[1], seg[2]; rtol=rtol)[1] - # int_n_e += seg_int - # # println("This segment: ", seg_int, "Total: ", int_n_e) - # end - # ch.n_e_line.data[ii] = int_n_e - integrand(s) = - n_e( - cp_n_e, - ep_n_e, - core_bnd, - SOLPS_bnd, - ep_space, - los(s)..., - ) .* dl_ds(s) - # function integrand!(y, s) - # n = Threads.nthreads() - # Threads.@threads for i ∈ 1:n - # y[i:n:end] .= integrand.(@view(s[i:n:end])) - # end - # end - # println("Calculating integrated electron density") - # ch.n_e_line.data[ii] = 0.0 - # for seg ∈ valid_segs - # @time ch.n_e_line.data[ii] += - # quadgk( - # BatchIntegrand{Float64}(integrand!), - # seg[1], - # seg[2], - # rtol, - # )[1] - # end ch.n_e_line.data[ii] = quadgk(integrand, 0, 1; rtol=rtol)[1] - # println(ch.n_e_line.data[ii]) - chord_in_core = get_intersections(cp_bnd, ep_space, chord_points, los) - core_chord_length = 0.0 - for seg ∈ chord_in_core - core_chord_length += - (seg[2] - seg[1]) * dl_ds((seg[1] + seg[2]) / 2) - end - # core_chord_integrand(s) = - # is_in_core( cp_bnd, ep_space, los(s)...) .* dl_ds(s) - # println("Calculating core chord length") - # @time core_chord_length = quadgk(core_chord_integrand, 0, 1, rtol)[1] - # println(core_chord_length) - ch.n_e_line_average.data[ii] = - ch.n_e_line.data[ii] / core_chord_length - # println("Average: ", ch.n_e_line_average.data[ii]) + 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] = @@ -290,6 +173,22 @@ function compute_interferometer(@nospecialize(ids::OMAS.dd), rtol::Float64=1e-1) end end +function get_sep_core_bnd(ep_grid_ggd) + ep_space = ep_grid_ggd.space[1] + core_bnd = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 15) + + core = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 22) + sol = SOLPS2IMAS.get_grid_subset_with_index(ep_grid_ggd, 23) + sep_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() + sep_bnd.element = + SOLPS2IMAS.subset_do( + intersect, + SOLPS2IMAS.get_subset_boundary(ep_space, sol), + SOLPS2IMAS.get_subset_boundary(ep_space, core), + ) + return sep_bnd, core_bnd +end + function rzphi2xyz( point::Union{OMAS.interferometer__channel___line_of_sight__first_point, OMAS.interferometer__channel___line_of_sight__second_point, @@ -304,25 +203,33 @@ function xyz2rz(x::Float64, y::Float64, z::Float64) return r, z end -function n_e( - cp_n_e::Function, - ep_n_e::Function, - core_bnd::OMAS.edge_profiles__grid_ggd___grid_subset, - SOLPS_bnd::OMAS.edge_profiles__grid_ggd___grid_subset, - ep_space::OMAS.edge_profiles__grid_ggd___space, - r::Float64, - z::Float64, -) - if (r, z) ∈ (core_bnd, ep_space) - return cp_n_e(r, z) - else# if (r, z) ∈ (SOLPS_bnd, ep_space) - return ep_n_e(r, z) - # else - # return 0.0 +function get_line_of_sight(points) + if length(points) == 2 + fp, sp = points + los = (s::Float64) -> xyz2rz((fp .+ s .* (sp .- fp))...) + dl_ds = (s::Float64) -> sqrt(sum((sp .- fp) .^ 2)) + else + fp, sp, tp = points + los = + (s::Float64) -> + if s <= 0.5 + return xyz2rz((fp .+ 2 * s .* (sp .- fp))...) + else + return xyz2rz((sp .+ 2 * (s - 0.5) .* (tp .- sp))...) + end + dl_ds = + (s::Float64) -> + if s <= 0.5 + return sqrt(sum((2 .* (sp .- fp)) .^ 2)) + else + return sqrt(sum((2 .* (sp .- tp)) .^ 2)) + end end + return los, dl_ds end -function get_intersections(subset, space, points, los) +function get_intersections(subset, space, points) + los, dl_ds = get_line_of_sight(points) nodes = space.objects_per_dimension[1].object edges = space.objects_per_dimension[2].object if length(points) == 2 @@ -355,13 +262,8 @@ function get_intersections(subset, space, points, los) filt_segs = Tuple{Float64, Float64}[] for ii ∈ 1:length(segs)-1 check_at = (segs[ii] + segs[ii+1]) / 2 - # print("Seg : ", los(segs[ii]), " to ", los(segs[ii+1]), ", Mid: ") - # print(check_at, ": ", los(check_at), ": ") if los(check_at) ∈ (subset, space) append!(filt_segs, [(segs[ii], segs[ii+1])]) - # println("Present") - # else - # println("Absent") end end return filt_segs @@ -397,3 +299,14 @@ function intersection_s( ) return intersection_s(fp, sp, xyz2rz(l1...), xyz2rz(l2...)) end + +function get_core_chord_length(sep_bnd, ep_space, chord_points) + los, dl_ds = get_line_of_sight(chord_points) + chord_in_core = get_intersections(sep_bnd, ep_space, chord_points) + core_chord_length = 0.0 + for seg ∈ chord_in_core + center_of_seg = (seg[1] + seg[2]) / 2 + core_chord_length += (seg[2] - seg[1]) * dl_ds(center_of_seg) + end + return core_chord_length +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 0baf7af..e52606a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,7 @@ using SynthDiag: add_interferometer! using SD4SOLPS: preparation using Test +using Printf @testset "interferometer" begin eqdsk_file = "g002296.00200" @@ -9,5 +10,23 @@ using Test ids = preparation(eqdsk_file, sample_dir; filename=test_dir * "/output") add_interferometer!("$(@__DIR__)/../src/default_interferometer.json", ids) # Just checking if the function runs through for now + for ch ∈ ids.interferometer.channel + println() + println("-"^49) + println("-"^49) + println("Channel $(ch.name)") + println("-"^49) + @printf("|%15s|%15s|%15s|\n", "time", "int_n_e", "n_e_average") + println("-"^49) + for ii ∈ eachindex(ch.n_e_line.data) + @printf( + "|%15.3e|%15.3e|%15.3e|\n", + ch.n_e_line.time[ii], + ch.n_e_line.data[ii], + ch.n_e_line_average.data[ii] + ) + end + println("-"^49) + end @test true end From 7df0790f9beb91d3d27f9f6d27c4cd8ff6a6a526 Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Mon, 9 Oct 2023 18:21:48 -0700 Subject: [PATCH 12/12] Added newline character at end of file --- src/interferometer.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/interferometer.jl b/src/interferometer.jl index 9c33833..e9dfe9e 100644 --- a/src/interferometer.jl +++ b/src/interferometer.jl @@ -309,4 +309,4 @@ function get_core_chord_length(sep_bnd, ep_space, chord_points) core_chord_length += (seg[2] - seg[1]) * dl_ds(center_of_seg) end return core_chord_length -end \ No newline at end of file +end