diff --git a/.gitignore b/.gitignore index 6b72f03..d31be6e 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ # committed for packages, but should be committed for applications that require a static # environment. Manifest.toml +examples/Project.toml diff --git a/Project.toml b/Project.toml index 2db44d5..073e5ca 100644 --- a/Project.toml +++ b/Project.toml @@ -10,7 +10,6 @@ Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" OMAS = "91cfaa06-6526-4804-8666-b540b3feef2f" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" -SOLPS2IMAS = "09becab6-0636-4c23-a92a-2b3723265c31" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/examples/Project.toml b/examples/Project.toml deleted file mode 100644 index 1cd060b..0000000 --- a/examples/Project.toml +++ /dev/null @@ -1,11 +0,0 @@ -name = "GGDUtils_examples" -authors = ["Anchal Gupta "] -version = "0.1.0" - -[deps] -GGDUtils = "b7b5e640-9b39-4803-84eb-376048795def" -GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" -LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" -OMAS = "91cfaa06-6526-4804-8666-b540b3feef2f" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -SOLPS2IMAS = "09becab6-0636-4c23-a92a-2b3723265c31" diff --git a/examples/plotting.ipynb b/examples/plotting.ipynb index 3ebd5f6..eabcd48 100644 --- a/examples/plotting.ipynb +++ b/examples/plotting.ipynb @@ -25,10 +25,7 @@ "using Pkg\n", "Pkg.activate(\"./\")\n", "Pkg.add(url=\"git@github.com:ProjectTorreyPines/OMAS.jl.git\")\n", - "Pkg.add(url=\"git@github.com:ProjectTorreyPines/SOLPS2IMAS.jl.git\", rev=\"dev\")\n", - "Pkg.add(path=\"../\")\n", - "Pkg.add(PackageSpec(name=\"GR\", version=\"0.72.9\"))\n", - "Pkg.pin(\"GR\")\n", + "Pkg.develop(path=\"../\")\n", "Pkg.add(\"Plots\")\n", "Pkg.add(\"LaTeXStrings\")" ] @@ -39,7 +36,7 @@ "metadata": {}, "outputs": [], "source": [ - "using SOLPS2IMAS\n", + "using OMAS\n", "using GGDUtils\n", "using Plots\n", "using LaTeXStrings" @@ -51,14 +48,9 @@ "metadata": {}, "outputs": [], "source": [ - "b2gmtry = \"../samples/b2fgmtry\"\n", - "b2output = \"../samples/b2time_red.nc\"\n", - "gsdesc = \"../samples/gridspacedesc.yml\"\n", - "b2mn = \"../samples/b2mn.dat\"\n", - "dd = solps2imas(b2gmtry, b2output, gsdesc, b2mn)\n", + "dd = OMAS.json2imas(\"../samples/time_dep_edge_profiles_last_step_only.json\")\n", "grid_ggd = dd.edge_profiles.grid_ggd[1]\n", - "space = grid_ggd.space[1]\n", - "dd.edge_profiles.ggd[1].electrons.density[5].grid_subset_index = 5 # This is a bug in SOLPS2IMAS to be fixed" + "space = grid_ggd.space[1]" ] }, { @@ -82,21 +74,21 @@ "\n", "# You can overlay any subset by giving a second argument\n", "# Labels \n", - "plot!(space, SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 6), markercolor=:chocolate1)\n", - "plot!(space, SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 7), linecolor=:red, linewidth=2)\n", - "plot!(space, SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 8), linecolor=:darkred, linewidth=2)\n", - "plot!(space, SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 9), linecolor=:limegreen, linewidth=2)\n", - "plot!(space, SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 10), linecolor=:darkgreen, linewidth=2)\n", - "plot!(space, SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 11), linecolor=:cyan, linewidth=2)\n", - "# plot!(space, SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 12), linecolor=:teal, linewidth=1)\n", - "plot!(space, SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 13), linecolor=:royalblue1, linewidth=2)\n", - "plot!(space, SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 14), linecolor=:navyblue, linewidth=2)\n", - "plot!(space, SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 15), linecolor=:fuchsia, linewidth=2, linestyle=:dash)\n", - "plot!(space, SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 16), linecolor=:purple4, linewidth=2, linestyle=:dash)\n", - "plot!(space, SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 101), markershape=:rect, markercolor=:royalblue1)\n", - "# plot!(space, SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 102), markershape=:rect, markercolor=:maroon)\n", - "plot!(space, SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 103), markershape=:diamond, markercolor=:fuchsia)\n", - "plot!(space, SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 104), markershape=:diamond, markercolor=:purple4)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 6), markercolor=:chocolate1)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 7), linecolor=:red, linewidth=2)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 8), linecolor=:darkred, linewidth=2)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 9), linecolor=:limegreen, linewidth=2)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 10), linecolor=:darkgreen, linewidth=2)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 11), linecolor=:cyan, linewidth=2)\n", + "# plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 12), linecolor=:teal, linewidth=1)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 13), linecolor=:royalblue1, linewidth=2)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 14), linecolor=:navyblue, linewidth=2)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 15), linecolor=:fuchsia, linewidth=2, linestyle=:dash)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 16), linecolor=:purple4, linewidth=2, linestyle=:dash)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 101), markershape=:rect, markercolor=:royalblue1)\n", + "# plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 102), markershape=:rect, markercolor=:maroon)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 103), markershape=:diamond, markercolor=:fuchsia)\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 104), markershape=:diamond, markercolor=:purple4)\n", "\n", "# Legend is supressed unless asked for specifically\n", "plot!(legend=true)\n", @@ -120,7 +112,8 @@ "gr() # Fast and can save pdf\n", "# plotlyjs() # Use for interactive plot, can only save png\n", "\n", - "plot(dd.edge_profiles.grid_ggd, dd.edge_profiles.ggd[1].electrons.density[5], colorbar_title=\"Electrons density / \" * L\"m^{-3}\")" + "n_e = GGDUtils.get_prop_with_grid_subset_index(dd.edge_profiles.ggd[1].electrons.density, 5)\n", + "plot(dd.edge_profiles.grid_ggd, n_e, colorbar_title=\"Electrons density / \" * L\"m^{-3}\")" ] }, { @@ -140,16 +133,9 @@ "gr() # Fast and can save pdf\n", "# plotlyjs() # Use for interactive plot, can only save png\n", "\n", - "plot(dd.edge_profiles.grid_ggd, dd.edge_profiles.ggd[1].electrons.density[5]) # Note default label in colorbar\n", - "plot!(space, SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 16), linecolor=:red, linewidth=2, linestyle=:solid, label=\"Separatix\", legend=true)" + "plot(dd.edge_profiles.grid_ggd, n_e) # Note default label in colorbar\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 16), linecolor=:red, linewidth=2, linestyle=:solid, label=\"Separatix\", legend=true)" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/examples/plotting_interferometer.ipynb b/examples/plotting_interferometer.ipynb new file mode 100644 index 0000000..917d0a5 --- /dev/null +++ b/examples/plotting_interferometer.ipynb @@ -0,0 +1,222 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Plotting interferometer using GGDUtils\n", + " \n", + " For running this notebook, you need to install package IJulia in your home environment (that is messy, but that is the only way I know right now). So in your terminal:\n", + " ```\n", + " % julia\n", + " julia > ]\n", + " (@v1.9 pkg) pkg> add IJulia\n", + " ```\n", + "\n", + " After this, Julia kernel would appear in your jupyter notebooks as an option. This also works for julia notebooks directly opened on VSCode. Select the Julia kernel to run this notebook." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using Pkg\n", + "Pkg.activate(\"./\")\n", + "Pkg.add(url=\"git@github.com:ProjectTorreyPines/OMAS.jl.git\")\n", + "Pkg.develop(path=\"../\")\n", + "Pkg.add(\"Plots\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using OMAS\n", + "using GGDUtils\n", + "using Plots" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Load the interferometer" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ids = OMAS.json2imas(\"$(@__DIR__)/../samples/time_dep_edge_profiles_with_interferometer.json\")\n", + "grid_ggd = ids.edge_profiles.grid_ggd[1]\n", + "space = grid_ggd.space[1]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plotting the interferometer geometry on top of SOLPS mesh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(space)\n", + "plot!(ids.interferometer) # Default plot_type is :los \n", + "plot!(legend=true)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can provide custom lengthand thickness of mirror to be plotted and linewidth of the laser beams" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(space)\n", + "plot!(ids.interferometer, mirror_length=0.7, linewidth=4, mirror_thickness=0.2)\n", + "plot!(legend=true)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Or you can choose to omit the mirror" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(space)\n", + "plot!(ids.interferometer, mirror=false)\n", + "plot!(legend=true)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can plot a single channel as well. You can override the in-built channel name for the label." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(space)\n", + "plot!(ids.interferometer.channel[1], label=\"Channel 1\")\n", + "plot!(legend=true)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plotting interferometer data vs time\n", + "\n", + " * Use plot_type=:n_e for integrated electron density data\n", + " * Use plot_type=:n_e_average for averaged electron density data\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(ids.interferometer, plot_type=:n_e)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(ids.interferometer, plot_type=:n_e_average)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Again, to plot an individual channel, just provide the channel with correct plot_type" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(ids.interferometer.channel[1], plot_type=:n_e_average)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.9.2", + "language": "julia", + "name": "julia-1.9" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.9.2" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/samples/.gitignore b/samples/.gitignore index 023a568..d203298 100644 --- a/samples/.gitignore +++ b/samples/.gitignore @@ -1,4 +1,2 @@ -/edge_profiles.h5 -/b2time.nc -/b2mn.dat -/b2fgmtry +/time_dep_edge_profiles_last_step_only.json +/time_dep_edge_profiles_with_interferometer.json diff --git a/samples/b2fgmtry.dvc b/samples/b2fgmtry.dvc deleted file mode 100644 index 0701492..0000000 --- a/samples/b2fgmtry.dvc +++ /dev/null @@ -1,12 +0,0 @@ -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/samples/b2mn.dat.dvc b/samples/b2mn.dat.dvc deleted file mode 100644 index e13ac93..0000000 --- a/samples/b2mn.dat.dvc +++ /dev/null @@ -1,12 +0,0 @@ -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/samples/b2time.nc.dvc b/samples/b2time.nc.dvc deleted file mode 100644 index b4a7624..0000000 --- a/samples/b2time.nc.dvc +++ /dev/null @@ -1,13 +0,0 @@ -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/samples/b2time_red.nc b/samples/b2time_red.nc deleted file mode 100644 index 91cd47f..0000000 Binary files a/samples/b2time_red.nc and /dev/null differ diff --git a/samples/edge_profiles.h5.dvc b/samples/edge_profiles.h5.dvc deleted file mode 100644 index 4af263a..0000000 --- a/samples/edge_profiles.h5.dvc +++ /dev/null @@ -1,12 +0,0 @@ -md5: 583339f63b0b4d2aa3f7f7a659a5b07c -frozen: true -deps: -- path: ITER_Lore_2296_00000/IMAS/edge_profiles.h5 - repo: - url: git@github.com:ProjectTorreyPines/SOLPSTestSamples.git - rev_lock: df499f1275428ec06175c48dc4af6ecbd5ec6117 -outs: -- md5: c3d60ce69e9d3a9b9ec621e9276858e0 - size: 168176600 - hash: md5 - path: edge_profiles.h5 diff --git a/samples/gridspacedesc.yml b/samples/gridspacedesc.yml deleted file mode 100644 index 0c82f2a..0000000 --- a/samples/gridspacedesc.yml +++ /dev/null @@ -1,196 +0,0 @@ -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/samples/time_dep_edge_profiles_last_step_only.json.dvc b/samples/time_dep_edge_profiles_last_step_only.json.dvc new file mode 100644 index 0000000..67f20fd --- /dev/null +++ b/samples/time_dep_edge_profiles_last_step_only.json.dvc @@ -0,0 +1,12 @@ +md5: 3ea764e4095ca2b8a5359cd76259e26a +frozen: true +deps: +- path: ITER_Lore_2296_00000/IMAS/time_dep_edge_profiles_last_step_only.json + repo: + url: git@github.com:ProjectTorreyPines/SOLPSTestSamples.git + rev_lock: 88f14955a52cd4d4fd7e86e2274436a34f4203ec +outs: +- md5: db47dee590648770a5ce49c0e0706acf + size: 4087084 + hash: md5 + path: time_dep_edge_profiles_last_step_only.json diff --git a/samples/time_dep_edge_profiles_with_interferometer.json.dvc b/samples/time_dep_edge_profiles_with_interferometer.json.dvc new file mode 100644 index 0000000..c45c30b --- /dev/null +++ b/samples/time_dep_edge_profiles_with_interferometer.json.dvc @@ -0,0 +1,12 @@ +md5: fb9792d06bb5d13f3a8803654b78624c +frozen: true +deps: +- path: ITER_Lore_2296_00000/IMAS/time_dep_edge_profiles_with_interferometer.json + repo: + url: git@github.com:ProjectTorreyPines/SOLPSTestSamples.git + rev_lock: 80688afd131b17455c3320a4485bbf2b7a55c80c +outs: +- md5: 610272ceabe4a62dd776e0f3c08d2036 + size: 40662828 + hash: md5 + path: time_dep_edge_profiles_with_interferometer.json diff --git a/src/GGDUtils.jl b/src/GGDUtils.jl index 67c639c..6c7fb25 100644 --- a/src/GGDUtils.jl +++ b/src/GGDUtils.jl @@ -2,6 +2,8 @@ module GGDUtils using OMAS: OMAS +const inv_16pi = 1.0 / (16π) + export project_prop_on_subset! export get_subset_centers export get_prop_with_grid_subset_index diff --git a/src/interpolations.jl b/src/interpolations.jl index f8992ef..1d4f9c2 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -68,11 +68,9 @@ minimizing bending energy of a surface. http://www.geometrictools.com/Documentation/ThinPlateSplines.pdf Eq(28) """ function _G(x1::Tuple{U, U}, x2::Tuple{U, U}) where {U <: Real} - r = sqrt(sum((x1 .- x2) .^ 2)) - if r == 0 - return 0 - end - return r^2 * log(r) / 8 / π + r2 = sum((x1 .- x2) .^ 2) # r² + # Note this uses log(r) / 8π = log(r²) / 16π + return (r2 == 0) ? zero(U) : r2 * log(r2) * inv_16pi end """ @@ -87,21 +85,15 @@ minimum absolute value. This is done to avoid numerical issues with interpolatio Return values are conditioned y and inverse conditioning function. """ function _condition_y(y::Vector{T}) where {T <: Real} - do_log = false ylims = extrema(y) - if prod(ylims) > 0 - if ylims[2] / ylims[1] > 100 - do_log = true - end - end - if do_log - norm_by = minimum(abs.(ylims)) * sign(ylims[1]) - return log10.(y ./ norm_by), (cy) -> (10 .^ (cy)) * norm_by - else - norm_by = ylims[2] - ylims[1] - mean_y = mean(y) - return (y .- mean_y) ./ norm_by, (cy) -> (cy .* norm_by) .+ mean_y + do_log = (prod(ylims) > 0) && (ylims[2] / ylims[1] > 100) + norm_by = do_log ? minimum(abs.(ylims)) * sign(ylims[1]) : ylims[2] - ylims[1] + mean_y = mean(y) + cy = do_log ? log10.(y ./ norm_by) : (y .- mean_y) ./ norm_by + inv_cy = let norm_by = norm_by, do_log = do_log, mean_y = mean_y + x -> do_log ? (10.0^x) * norm_by : (x * norm_by) + mean_y end + return cy, inv_cy end """ @@ -129,6 +121,12 @@ function get_TPS_mats(x::Vector{Tuple{U, U}}) where {U <: Real} return Minv, N, (N' * Minv * N)^(-1) * N' * Minv, x end +function get_interp_val(r, z, x, a, b, inv_cy) + tot = sum(a[k] * _G((r, z), x[k]) for k ∈ eachindex(a)) + tot += b[1] + r * b[2] + z * b[3] + return inv_cy(tot) +end + """ interp( y::Vector{T}, @@ -147,11 +145,10 @@ function interp( # From Eq(31) b = y2b * cy a = Minv * (cy - N * b) - function get_interp_val(r::Real, z::Real) - return inv_cy(sum(a .* [_G((r, z), xi) for xi ∈ x]) + sum(b .* [1, r, z])) + return let x = x, a = a, b = b, inv_cy = inv_cy + g(r, z) = get_interp_val(r, z, x, a, b, inv_cy) + g(gp::Tuple{V, V}) where {V <: Real} = g(gp...) end - get_interp_val(gp::Tuple{V, V}) where {V <: Real} = get_interp_val(gp...) - return get_interp_val end """ @@ -322,10 +319,11 @@ function interp( prop_arr::Vector{T}, TPS_mats::Tuple{Matrix{U}, Matrix{U}, Matrix{U}, Vector{Tuple{U, U}}}, grid_subset_index::Int, - value_field::Symbol=:values, -) where {T <: edge_profiles__prop_on_subset, U <: Real} + value_field::Val{V}=Val(:values), +) where {T <: edge_profiles__prop_on_subset, U <: Real, V} prop = get_prop_with_grid_subset_index(prop_arr, grid_subset_index) - return interp(getfield(prop, value_field), TPS_mats) + field = getfield(prop, V) + return interp(field, TPS_mats) end const RHO_EXT_POS = [1.0001, 1.1, 5] @@ -363,9 +361,10 @@ function interp(eqt::OMAS.equilibrium__time_slice) prepend!(rhon_eq_ext, RHO_EXT_NEG) rz2psin = linear_interpolation((r_eq, z_eq), psinrz) psin2rhon = linear_interpolation(psin_eq_ext, rhon_eq_ext) - get_interp_val(r::Real, z::Real) = psin2rhon(rz2psin(r, z)) - get_interp_val(rz::Tuple{Real, Real}) = get_interp_val(rz...) - return get_interp_val + g = let psin2rhon = psin2rhon, rz2psin = rz2psin + (r, z) -> psin2rhon(rz2psin(r, z)) + end + return g end """ @@ -423,9 +422,10 @@ function interp( rz2rho::Function, ) where {T <: Real} itp = interp(prop, prof) - get_interp_val(r::Real, z::Real) = itp.(rz2rho(r, z)) - get_interp_val(rz::Tuple{Real, Real}) = get_interp_val(rz...) - return get_interp_val + g = let itp = itp + (r, z) -> itp(rz2rho(r, z)) + end + return g end """ diff --git a/src/recipes.jl b/src/recipes.jl index a2bacce..3e4ef79 100644 --- a/src/recipes.jl +++ b/src/recipes.jl @@ -1,5 +1,6 @@ using RecipesBase using ColorSchemes: ColorSchemes +import Statistics: norm, dot @recipe function f(space::OMAS.edge_profiles__grid_ggd___space) nodes = space.objects_per_dimension[1].object @@ -46,15 +47,15 @@ end yaxis --> "Z / m" subset_edge_inds = [] label_assigned = false - if subset.element[1].object[1].dimension == 2 + if subset.element[1].object[1].dimension == 3 for ele ∈ subset.element for bnd_ind ∈ cells[ele.object[1].index].boundary union!(subset_edge_inds, bnd_ind) end end - elseif subset.element[1].object[1].dimension == 1 + elseif subset.element[1].object[1].dimension == 2 subset_edge_inds = [ele.object[1].index for ele ∈ subset.element] - elseif subset.element[1].object[1].dimension == 0 + elseif subset.element[1].object[1].dimension == 1 for ele ∈ subset.element @series begin seriestype := :scatter @@ -100,14 +101,13 @@ end subset = get_grid_subset_with_index(grid_ggd, prop.grid_subset_index) space = grid_ggd.space[subset.element[1].object[1].space] nodes = space.objects_per_dimension[1].object - # edges = space.objects_per_dimension[2].object cells = space.objects_per_dimension[3].object legend --> false size --> [600, 900] xaxis --> "R / m" yaxis --> "Z / m" layout := @layout [a{0.95w} b] - if subset.element[1].object[1].dimension == 2 + if subset.element[1].object[1].dimension == 3 if :seriescolor in keys(plotattributes) color_scheme = plotattributes[:seriescolor] else @@ -176,3 +176,151 @@ end ) end end + +@recipe f( + ifo::OMAS.interferometer, +) = + for ch ∈ ifo.channel + @series begin + return ch + end + end + +@recipe function f( + ifo_ch::OMAS.interferometer__channel, +) + if :plot_type ∈ keys(plotattributes) + plot_type = plotattributes[:plot_type] + else + plot_type = :los + end + if plot_type == :los + @series begin + label --> ifo_ch.name + ifo_ch.line_of_sight + end + elseif plot_type == :n_e || plot_type == :n_e_line + @series begin + label --> ifo_ch.name + ifo_ch.n_e_line + end + elseif plot_type == :n_e_average || plot_type == :n_e_line_average + @series begin + label --> ifo_ch.name + ifo_ch.n_e_line_average + end + else + error("Invalid plot type. Choose between los (line of sight, R-Z plane), ", + "n_e (integrated n_e along line of sight vs time), ", + "and n_e_average (average n_e vs time)") + end +end + +@recipe function f( + ifo_ch_los::OMAS.interferometer__channel___line_of_sight, +) + size --> [600, 900] + xaxis --> "R / m" + yaxis --> "Z / m" + fp = ifo_ch_los.first_point + sp = ifo_ch_los.second_point + tp = ifo_ch_los.third_point + if tp == OMAS.interferometer__channel___line_of_sight__third_point() + tp = fp + end + + if :mirror ∈ keys(plotattributes) + mirror = plotattributes[:mirror] + else + mirror = true + end + if mirror + # Calculate mirror points + fsuv = Array([sp.r - fp.r, sp.z - fp.z]) + stuv = Array([tp.r - sp.r, tp.z - sp.z]) + fsuv = fsuv / norm(fsuv) + stuv = stuv / norm(stuv) + mirror_uv = fsuv + stuv + if mirror_uv[1] == 0 && mirror_uv[2] == 0 + mirror_uv = [fsuv[2], -fsuv[1]] + end + mirror_uv = mirror_uv / norm(mirror_uv) + if :mirror_length ∈ keys(plotattributes) + mirror_length = plotattributes[:mirror_length] + else + mirror_length = 0.5 + end + mirror_perp = ([-mirror_uv[2], mirror_uv[1]]) + mirror_perp = mirror_perp * sign(dot(mirror_perp, fsuv)) + if :mirror_thickness ∈ keys(plotattributes) + mirror_thickness = plotattributes[:mirror_thickness] + else + mirror_thickness = 0.1 + end + mirror_r1 = sp.r - mirror_uv[1] * mirror_length / 2 + mirror_r2 = sp.r + mirror_uv[1] * mirror_length / 2 + mirror_r3 = mirror_r2 + mirror_perp[1] * mirror_thickness + mirror_r4 = mirror_r1 + mirror_perp[1] * mirror_thickness + mirror_z1 = sp.z - mirror_uv[2] * mirror_length / 2 + mirror_z2 = sp.z + mirror_uv[2] * mirror_length / 2 + mirror_z3 = mirror_z2 + mirror_perp[2] * mirror_thickness + mirror_z4 = mirror_z1 + mirror_perp[2] * mirror_thickness + end + + # Draw line of sight + @series begin + seriestype := :path + linewidth --> 2 + [fp.r, sp.r, tp.r], [fp.z, sp.z, tp.z] + end + + if mirror + # Draw mirror + @series begin + seriestype := :path + label := "" + linecolor := :gray + linewidth := 2 + [mirror_r1, mirror_r2], [mirror_z1, mirror_z2] + end + @series begin + seriestype := :shape + label := "" + linealpha := 0 + fillstyle := :/ + fillcolor := :black + [mirror_r1, mirror_r2, mirror_r3, mirror_r4], + [mirror_z1, mirror_z2, mirror_z3, mirror_z4] + end + end +end + +@recipe function f( + ifo_ch_n_e_line::OMAS.interferometer__channel___n_e_line, +) + if :average ∈ keys(plotattributes) + @series begin + ifo_ch.n_e_line_average + end + else + xaxis --> "time / s" + yaxis --> "Integrrated n_e / m^-2" + @series begin + seriestype := :path + linewidth --> 2 + ifo_ch_n_e_line.time, ifo_ch_n_e_line.data + end + end +end + +@recipe function f( + ifo_ch_n_e_line_average::OMAS.interferometer__channel___n_e_line_average, +) + xaxis --> "time / s" + yaxis --> "Average n_e / m^-3" + @series begin + seriestype := :path + linewidth --> 2 + ifo_ch_n_e_line_average.time, ifo_ch_n_e_line_average.data + end +end diff --git a/src/subset_tools.jl b/src/subset_tools.jl index 5b269dd..1638574 100644 --- a/src/subset_tools.jl +++ b/src/subset_tools.jl @@ -1,4 +1,197 @@ -import SOLPS2IMAS: get_subset_space, get_grid_subset_with_index, get_subset_boundary +""" + add_subset_element!( + subset::OMAS.edge_profiles__grid_ggd___grid_subset, + sn::Int, + dim::Int, + index::Int, + in_subset=(x...) -> true; + kwargs..., + +) + +Adds a new element to gird_subset with properties space number (sn), dimension (dim), +and index (index). The element is added only if the function in_subset returns true. +""" +function add_subset_element!( + subset::OMAS.edge_profiles__grid_ggd___grid_subset, + sn::Int, + dim::Int, + index::Int, + in_subset=(x...) -> true; + kwargs..., +) + if in_subset(; kwargs...) + dd_ind = length(subset.element) + 1 + resize!(subset.element, dd_ind) + resize!(subset.element[dd_ind].object, 1) + subset.element[dd_ind].object[1].space = sn + subset.element[dd_ind].object[1].dimension = dim + subset.element[dd_ind].object[1].index = index + end +end + +""" + add_subset_element!( + subset, + sn, + dim, + index::Vector{Int}, + in_subset=(x...) -> true; + kwargs..., + +) + +Overloaded to work differently (faster) with list of indices to be added. +""" +function add_subset_element!( + subset::OMAS.edge_profiles__grid_ggd___grid_subset, + sn::Int, + dim::Int, + index::Vector{Int}, + in_subset=(x...) -> true; + kwargs..., +) + if in_subset(; kwargs...) + dd_start_ind = length(subset.element) + 1 + resize!(subset.element, length(subset.element) + length(index)) + dd_stop_ind = length(subset.element) + for (ii, dd_ind) ∈ enumerate(dd_start_ind:dd_stop_ind) + resize!(subset.element[dd_ind].object, 1) + subset.element[dd_ind].object[1].space = sn + subset.element[dd_ind].object[1].dimension = dim + subset.element[dd_ind].object[1].index = index[ii] + end + end +end + +""" + get_subset_space(space::OMAS.edge_profiles__grid_ggd___space, + elements::Vector{OMAS.edge_profiles__grid_ggd___grid_subset___element}) + +Returns an array of space object indices corresponding to the correct +objects_per_dimension (nodes, edges or cells) for the subset elements. +""" +function get_subset_space(space::OMAS.edge_profiles__grid_ggd___space, + elements::Vector{OMAS.edge_profiles__grid_ggd___grid_subset___element}) + nD = elements[1].object[1].dimension + nD_objects = space.objects_per_dimension[nD].object + return [nD_objects[ele.object[1].index] for ele ∈ elements] +end + +""" + get_grid_subset_with_index( + grid_ggd::OMAS.edge_profiles__grid_ggd, + grid_subset_index::Int, + +) + +Returns the grid_subset in a grid_ggd with the matching grid_subset_index +""" +function get_grid_subset_with_index( + grid_ggd::OMAS.edge_profiles__grid_ggd, + grid_subset_index::Int, +) + for subset ∈ grid_ggd.grid_subset + if subset.identifier.index == grid_subset_index + return subset + end + end + # BCL 12/8: Creates type instability, but maybe okay since it's a "simple" Union + # subset::Union{Int64, OMAS.edge_profiles__grid_ggd___grid_subset} + # + # Better would be to immediately throw an error or return nothing + return 0 # Indicates failure +end + +""" + get_subset_boundary_inds( + space::OMAS.edge_profiles__grid_ggd___space, + subset::OMAS.edge_profiles__grid_ggd___grid_subset, + +) + +Returns an array of space object indices corresponding to the boundary of the subset. +That means, it returns indices of nodes that are at the end of open edge subset or +it returns the indices of edges that are the the boundary of a cell subset. +Returns an empty array if the subset is 1D (nodes). +""" +function get_subset_boundary_inds( + space::OMAS.edge_profiles__grid_ggd___space, + subset::OMAS.edge_profiles__grid_ggd___grid_subset, +) + nD = subset.element[1].object[1].dimension + if nD > 1 # Only 2D (edges) and 3D (cells) subsets have boundaries + nD_objects = space.objects_per_dimension[nD].object + elements = [nD_objects[ele.object[1].index] for ele ∈ subset.element] + boundary_inds = Int[] + for ele ∈ elements + symdiff!(boundary_inds, [bnd.index for bnd ∈ ele.boundary]) + end + return boundary_inds + end + return [] # 1D (nodes) subsets have no boundary +end + +""" + get_subset_boundary( + space::OMAS.edge_profiles__grid_ggd___space, + subset::OMAS.edge_profiles__grid_ggd___grid_subset, + +) + +Returns an array of elements of grid_subset generated from the boundary of the subset +provided. The dimension of these elments is reduced by 1. +""" +function get_subset_boundary( + space::OMAS.edge_profiles__grid_ggd___space, + subset::OMAS.edge_profiles__grid_ggd___grid_subset, +) + ret_subset = OMAS.edge_profiles__grid_ggd___grid_subset() + boundary_inds = get_subset_boundary_inds(space, subset) + bnd_dim = subset.element[1].object[1].dimension - 1 + space_number = subset.element[1].object[1].space + add_subset_element!(ret_subset, space_number, bnd_dim, boundary_inds) + return ret_subset.element +end + +""" + subset_do(set_operator, + itrs::Vararg{Vector{OMAS.edge_profiles__grid_ggd___grid_subset___element}}; + space::OMAS.edge_profiles__grid_ggd___space=OMAS.edge_profiles__grid_ggd___space(), + use_nodes=false) + +Function to perform any set operation (intersect, union, setdiff etc.) on +subset.element to generate a list of elements to go to subset object. If use_nodes is +true, the set operation will be applied on the set of nodes from subset.element, space +argument is required for this. +Note: that the arguments are subset.element (not the subset itself). Similarly, the +return object is a list of OMAS.edge_profiles__grid_ggd___grid_subset___element. +""" +function subset_do(set_operator, + itrs::Vararg{Vector{OMAS.edge_profiles__grid_ggd___grid_subset___element}}; + space::OMAS.edge_profiles__grid_ggd___space=OMAS.edge_profiles__grid_ggd___space(), + use_nodes=false) + if use_nodes + ele_inds = set_operator( + [ + union([ + obj.nodes for obj ∈ get_subset_space(space, set_elements) + ]... + ) for set_elements ∈ itrs + ]..., + ) + dim = 1 + else + ele_inds = set_operator( + [[ele.object[1].index for ele ∈ set_elements] for set_elements ∈ itrs]..., + ) + dim = itrs[1][1].object[1].dimension + end + ret_subset = OMAS.edge_profiles__grid_ggd___grid_subset() + space_number = itrs[1][1].object[1].space + add_subset_element!(ret_subset, space_number, dim, ele_inds) + return ret_subset.element +end """ get_subset_centers(space::OMAS.edge_profiles__grid_ggd___space, @@ -8,11 +201,12 @@ Returns an array of tuples corresponding to (r,z) coordinates of the center of cells or the center of edges in the subset space. """ function get_subset_centers(space::OMAS.edge_profiles__grid_ggd___space, - subset::OMAS.edge_profiles__grid_ggd___grid_subset) + subset::OMAS.edge_profiles__grid_ggd___grid_subset, +) subset_space = get_subset_space(space, subset.element) grid_nodes = space.objects_per_dimension[1].object return [ - Tuple(mean([grid_nodes[node].geometry for node ∈ obj.nodes])) for + Tuple(mean(SVector{2}(grid_nodes[node].geometry) for node ∈ obj.nodes)) for obj ∈ subset_space ] end @@ -23,7 +217,7 @@ end from_subset::OMAS.edge_profiles__grid_ggd___grid_subset, to_subset::OMAS.edge_profiles__grid_ggd___grid_subset, space::OMAS.edge_profiles__grid_ggd___space, -) +) This function can be used to add another instance on a property vector representing the value in a new subset that can be taken as a projection from an existing larger subset. @@ -64,7 +258,7 @@ function project_prop_on_subset!(prop_arr::Vector{T}, to_subset::OMAS.edge_profiles__grid_ggd___grid_subset, space::OMAS.edge_profiles__grid_ggd___space, value_field::Symbol=:values; - interp_method=:thin_plate_spline, + interp_method::Symbol=:thin_plate_spline, interp_kwargs=Dict(), ) where {T <: edge_profiles__prop_on_subset} if from_subset.element[1].object[1].dimension == @@ -161,10 +355,10 @@ end Overloading ∈ operator to check if a point is inside a subset of space. -If the subset is 0-dimensional, all points are searched. If the subset is 1-dimensional, +If the subset is 1-dimensional, all points are searched. If the subset is 2-dimensional, it is checked if the point is within the enclosed area. It is assumed that a -1-dimensional subset used in such a context will form a closed area. If the subset is -2-dimensional, its boundary is calculated on the fly. If used multiple times, it is +2-dimensional subset used in such a context will form a closed area. If the subset is +3-dimensional, its boundary is calculated on the fly. If used multiple times, it is recommended to calculate the boundary once and store it in a variable. """ function Base.:∈( @@ -179,12 +373,12 @@ function Base.:∈( dim = subset.element[1].object[1].dimension nodes = space.objects_per_dimension[1].object edges = space.objects_per_dimension[2].object - if dim == 2 + if dim == 3 subset_bnd = OMAS.edge_profiles__grid_ggd___grid_subset() subset_bnd.element = get_subset_boundary(space, subset) - elseif dim == 1 + elseif dim == 2 subset_bnd = subset - elseif dim == 0 + elseif dim == 1 for ele ∈ subset.element node = nodes[ele.object[1].index] if node.geometry[1] == r && node.geometry[2] == z @@ -199,21 +393,17 @@ function Base.:∈( count = 0 for ele ∈ subset_bnd.element edge = edges[ele.object[1].index] - r_max = maximum([nodes[node].geometry[1] for node ∈ edge.nodes]) - r_min = minimum([nodes[node].geometry[1] for node ∈ edge.nodes]) + r_max = maximum(nodes[node].geometry[1] for node ∈ edge.nodes) + r_min = minimum(nodes[node].geometry[1] for node ∈ edge.nodes) if r_min <= r < r_max - z_max = maximum([nodes[node].geometry[2] for node ∈ edge.nodes]) + z_max = maximum(nodes[node].geometry[2] for node ∈ edge.nodes) if z < z_max count += 1 end end end # If it is even, the point is outside the boundary - if count % 2 == 1 - return true - else - return false - end + return count % 2 == 1 end function get_prop_with_grid_subset_index( diff --git a/test/runtests.jl b/test/runtests.jl index bc8d69a..c0c2b74 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,6 @@ -import GGDUtils: interp, get_kdtree, project_prop_on_subset! -import OMAS: h5i2imas +import GGDUtils: interp, get_kdtree, project_prop_on_subset!, get_grid_subset_with_index +import OMAS: json2imas import Statistics: mean -import SOLPS2IMAS: solps2imas, get_grid_subset_with_index using Test using ArgParse: ArgParse @@ -34,6 +33,11 @@ function parse_commandline() end args = parse_commandline() +print("json2imas() time: ") +@time ids = json2imas( + "$(@__DIR__)/../samples/time_dep_edge_profiles_last_step_only.json", +) + if args["interp"] @testset "interp" begin # ids = h5i2imas("$(@__DIR__)/../samples/edge_profiles.h5") @@ -41,7 +45,6 @@ if args["interp"] b2output = "$(@__DIR__)/../samples/b2time.nc" gsdesc = "$(@__DIR__)/../samples/gridspacedesc.yml" b2mn = "$(@__DIR__)/../samples/b2mn.dat" - ids = solps2imas(b2gmtry, b2output, gsdesc, b2mn) n_e = ids.edge_profiles.ggd[1].electrons.density[1] grid_ggd = ids.edge_profiles.grid_ggd[1] space = grid_ggd.space[1] @@ -54,7 +57,8 @@ if args["interp"] grid_val = ids.edge_profiles.ggd[1].electrons.density[1].values[chosen_index] # test interp(prop, grid_ggd) - get_n_e = interp(n_e, grid_ggd) + print("interp(prop, grid_ggd) time: ") + @time get_n_e = interp(n_e, grid_ggd) searched_val = get_n_e(cell_center...) println("Electron density at: ", cell_center) println("Grid Value: ", grid_val) @@ -63,21 +67,26 @@ if args["interp"] # test interp(prop_arr, space, subset) subset = get_grid_subset_with_index(grid_ggd, 5) - get_n_e = interp(ids.edge_profiles.ggd[1].electrons.density, space, subset) + print("interp(prop_arr, space, subset) time: ") + @time get_n_e = + interp(ids.edge_profiles.ggd[1].electrons.density, space, subset) searched_val = get_n_e(cell_center...) @test abs.((grid_val .- searched_val) ./ grid_val) < allowed_rtol # test interp(prop_arr, grid_ggd, grid_subset_index) - get_n_e = interp(ids.edge_profiles.ggd[1].electrons.density, grid_ggd, 5) + print("interp(prop_arr, grid_ggd, grid_subset_index) time: ") + @time get_n_e = interp(ids.edge_profiles.ggd[1].electrons.density, grid_ggd, 5) searched_val = get_n_e(cell_center...) @test abs.((grid_val .- searched_val) ./ grid_val) < allowed_rtol # Use the kdtree to interpolate several quantities using # inverse distance weighing kdtree = get_kdtree(space) - get_T_e = + print("interp(prop_arr, kdtree) time: ") + @time get_T_e = interp(ids.edge_profiles.ggd[1].electrons.temperature[1].values, kdtree) - get_n_e = interp(ids.edge_profiles.ggd[1].electrons.density[1].values, kdtree) + @time get_n_e = + interp(ids.edge_profiles.ggd[1].electrons.density[1].values, kdtree) grid_val = ids.edge_profiles.ggd[1].electrons.temperature[1].values[chosen_index] @@ -116,27 +125,24 @@ end if args["projection"] @testset "project_prop_on_subset!" begin - b2gmtry = "$(@__DIR__)/../samples/b2fgmtry" - b2output = "$(@__DIR__)/../samples/b2time.nc" - gsdesc = "$(@__DIR__)/../samples/gridspacedesc.yml" - b2mn = "$(@__DIR__)/../samples/b2mn.dat" - dd = solps2imas(b2gmtry, b2output, gsdesc, b2mn) - space = dd.edge_profiles.grid_ggd[1].space[1] - prop = dd.edge_profiles.ggd[1].electrons.density + space = ids.edge_profiles.grid_ggd[1].space[1] + prop = ids.edge_profiles.ggd[1].electrons.density # All cells - from_subset = get_grid_subset_with_index(dd.edge_profiles.grid_ggd[1], 5) - kdtree = get_kdtree(space, from_subset) + from_subset = get_grid_subset_with_index(ids.edge_profiles.grid_ggd[1], 5) # separatix - to_subset = get_grid_subset_with_index(dd.edge_profiles.grid_ggd[1], 16) - separatix_centers, values_at_separatix = + to_subset = get_grid_subset_with_index(ids.edge_profiles.grid_ggd[1], 16) + print("project_prop_on_subset!(prop, from_subset, to_subset, space) time: ") + @time separatix_centers, values_at_separatix = project_prop_on_subset!(prop, from_subset, to_subset, space) # println("Projected to separatix:") # for ii ∈ eachindex(separatix_centers) # println(separatix_centers[ii], ": ", values_at_separatix[ii]) # end - subset_core = get_grid_subset_with_index(dd.edge_profiles.grid_ggd[1], 22) - core_element_inds, values_at_core = + subset_core = + get_grid_subset_with_index(ids.edge_profiles.grid_ggd[1], 22) + print("project_prop_on_subset!(prop, from_subset, subset_core) time: ") + @time core_element_inds, values_at_core = project_prop_on_subset!(prop, from_subset, subset_core) # println("Project to core:") # for ii ∈ eachindex(core_element_inds) @@ -153,12 +159,7 @@ end if args["in"] @testset "test ∈" begin - b2gmtry = "$(@__DIR__)/../samples/b2fgmtry" - b2output = "$(@__DIR__)/../samples/b2time.nc" - gsdesc = "$(@__DIR__)/../samples/gridspacedesc.yml" - b2mn = "$(@__DIR__)/../samples/b2mn.dat" - dd = solps2imas(b2gmtry, b2output, gsdesc, b2mn) - grid_ggd = dd.edge_profiles.grid_ggd[1] + grid_ggd = ids.edge_profiles.grid_ggd[1] space = grid_ggd.space[1] subset_corebnd = get_grid_subset_with_index(grid_ggd, 15) subset_sol = get_grid_subset_with_index(grid_ggd, 23)