diff --git a/Project.toml b/Project.toml index 8b2bc92..a3e1890 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SD4SOLPS" uuid = "d8db6f1b-e564-4c04-bba3-ef399f78c070" authors = ["David Eldon "] -version = "0.1.0" +version = "0.2.0" [deps] ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" @@ -9,10 +9,12 @@ Contour = "d38c429a-6771-53c6-b99e-75d170b6e991" EFIT = "cda752c5-6b03-55a3-9e33-132a441b0c17" GGDUtils = "b7b5e640-9b39-4803-84eb-376048795def" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" +JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" OMAS = "91cfaa06-6526-4804-8666-b540b3feef2f" PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab" PlotUtils = "995b91a9-d308-5afd-9ec6-746e21dbc043" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +PolygonOps = "647866c9-e3ac-4575-94e7-e3d426903924" SOLPS2IMAS = "09becab6-0636-4c23-a92a-2b3723265c31" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" diff --git a/README.md b/README.md index b4143c7..ed84d70 100644 --- a/README.md +++ b/README.md @@ -9,3 +9,55 @@ Steps: 2) Load equilibrium (that the SOLPS mesh was based on) into IMAS DD format 3) Make assumptions to extend profiles into the core and far SOL, if needed 4) Run synthetic diagnostic models and record output + + +## Installation + +### Using make file + +#### Option 1: Download the [example directory](https://github.com/ProjectTorreyPines/SD4SOLPS.jl/tree/master/example) in this repo: + +```bash +% make help +Help Menu + +make env_with_cloned_repo: Creates a Julia environment with the cloned repositories +make env_with_git_url: Creates a Julia environment with the git urls without creating local clones +make clean: Deletes Project.toml and Manifest.toml for a fresh start +make Clean: Deletes Project.toml, Manifest.toml, and any cloned repositories for a fresh start +``` +##### Environment with cloned repos + +`make env_with_cloned_repo` will clone all required git repos one level up from your current working directory and use them to define the environment files that will be generated in your working directory. This way, in future you can do `git pull` or change branches to update your working environment. + +##### Static environment with git url + +`make env_with_git_url` will not create local git clones for you, but it will still clone master branches of all the required repos in julia package management directory and will create a static environment for you that will remain tied to the version of packages when this environment was created. + +##### Clean up + +`make clean` or `make Clean` can be used clean up environment files and redo the environment generation. + +#### Option 2: Clone this repo + +```bash +% git clone git@github.com:ProjectTorreyPines/SD4SOLPS.jl.git +% cd example +$ make help +Help Menu + +make env_with_cloned_repo: Creates a Julia environment with the cloned repositories +make env_with_git_url: Creates a Julia environment with the git urls without creating local clones +make clean: Deletes Project.toml and Manifest.toml for a fresh start +make Clean: Deletes Project.toml, Manifest.toml, and any cloned repositories for a fresh start +``` + +Further options are same as above except for the difference that in case of cloning local copies of repos, they will be kept on same level as where you cloned SD4SOLPS.jl repo. + +### Manual Install + +Refer to the instructions on this [wiki page](https://github.com/ProjectTorreyPines/SD4SOLPS.jl/wiki). + +## Use + +Refer to the instructions on this [wiki page](https://github.com/ProjectTorreyPines/SD4SOLPS.jl/wiki/Demo). diff --git a/apply_format.sh b/apply_format.sh new file mode 100755 index 0000000..0101490 --- /dev/null +++ b/apply_format.sh @@ -0,0 +1,2 @@ +#!/bin/bash +julia -e 'using JuliaFormatter; format(".", verbose=true)' diff --git a/data/.gitignore b/data/.gitignore new file mode 100644 index 0000000..f6bc585 --- /dev/null +++ b/data/.gitignore @@ -0,0 +1,2 @@ +*.mesh_ext.json + diff --git a/example/create_env_with_clone.sh b/example/create_env_with_clone.sh new file mode 100644 index 0000000..f1c4225 --- /dev/null +++ b/example/create_env_with_clone.sh @@ -0,0 +1,38 @@ +#!/bin/bash + +is_repo=$(git rev-parse --is-inside-work-tree 2>/dev/null) +if [[ ${is_repo} == "true" ]] +then + SD4SOLPS_path="$(git rev-parse --show-toplevel)" + project_path="$(dirname ${SD4SOLPS_path})" +else + SD4SOLPS_path="" + project_path="$(dirname $(pwd))" +fi +echo "Cloning git repositories to "$project_path +orig_path=$(pwd) +cd $project_path +if [[ -z $SD4SOLPS_path ]] +then + git clone git@github.com:ProjectTorreyPines/SD4SOLPS.jl.git + SD4SOLPS_path=$project_path"/SD4SOLPS.jl" +fi +# in a for loop clone SynthDiag.jl SOLPS2IMAS.jl OMAS.jl GGDUtils.jl +for repo in SynthDiag.jl SOLPS2IMAS.jl OMAS.jl GGDUtils.jl +do + if [[ ! -d $repo ]] + then + git clone "git@github.com:ProjectTorreyPines/"$repo".git" + fi +done +cd $orig_path + +echo "Generating Project.toml and Manifest.toml" +OMAS_path=$project_path"/OMAS.jl" +EFIT_url="git@github.com:JuliaFusion/EFIT.jl.git" +julia --project=. -e 'using Pkg; Pkg.develop(path="'$OMAS_path'"); Pkg.instantiate()' +julia --project=. -e 'using Pkg; Pkg.add(url="'$EFIT_url'", rev="master"); Pkg.instantiate()' +for repo in GGDUtils SOLPS2IMAS SynthDiag SD4SOLPS; do + repo_path=$project_path"/"$repo".jl" + julia --project=. -e 'using Pkg; Pkg.develop(path="'${repo_path}'"); Pkg.instantiate()' +done \ No newline at end of file diff --git a/example/create_env_with_git.sh b/example/create_env_with_git.sh new file mode 100644 index 0000000..38188f9 --- /dev/null +++ b/example/create_env_with_git.sh @@ -0,0 +1,11 @@ +#!/bin/bash + +echo "Generating Project.toml and Manifest.toml" +OMAS_url="git@github.com:ProjectTorreyPines/OMAS.jl.git" +EFIT_url="git@github.com:JuliaFusion/EFIT.jl.git" +julia --project=. -e 'using Pkg; Pkg.add(url="'$OMAS_url'", rev="master"); Pkg.instantiate()' +julia --project=. -e 'using Pkg; Pkg.add(url="'$EFIT_url'", rev="master"); Pkg.instantiate()' +for repo in GGDUtils SOLPS2IMAS SynthDiag SD4SOLPS; do + repo_url="git@github.com:ProjectTorreyPines/"$repo".jl.git" + julia --project=. -e 'using Pkg; Pkg.add(url="'$repo_url'", rev="master"); Pkg.instantiate()' +done \ No newline at end of file diff --git a/example/demo.ipynb b/example/demo.ipynb new file mode 100644 index 0000000..6841ccf --- /dev/null +++ b/example/demo.ipynb @@ -0,0 +1,461 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using Pkg;\n", + "Pkg.activate(\"./\")\n", + "Pkg.instantiate()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Loading SOLPS output into IMAS data structure" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using SOLPS2IMAS: SOLPS2IMAS" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "b2gmtry = \"samples/b2fgmtry\"\n", + "b2output = \"samples/b2time.nc\"\n", + "gsdesc = \"samples/gridspacedesc.yml\"\n", + "b2mn = \"samples/b2mn.dat\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ids = SOLPS2IMAS.solps2imas(b2gmtry, b2output, gsdesc, b2mn)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Visualising some properties" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using GGDUtils: GGDUtils\n", + "using Plots" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "IMAS ids store mesh information for edge profiles in grid_ggd. There are options to have multiple space representaions but typically you will have only one space describing the SOLPS mesh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "grid_ggd = ids.edge_profiles.grid_ggd[1] # First grid_ggd time slice. It is allowed to vary in time\n", + "space = grid_ggd.space[1] # First space in this grid_ggd" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plotting grid and subsets" + ] + }, + { + "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) # Simply plot the grid described in space, all common arguments to plot can be given here\n", + "\n", + "# You can overlay any subset by giving a second argument\n", + "# Labels \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", + "# Default labels are subset.identifier.name but can be changed by providing a label argument\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plotting 2D quantities as heatmaps" + ] + }, + { + "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", + "n_e = GGDUtils.get_prop_with_grid_subset_index(ids.edge_profiles.ggd[1].electrons.density, 5)\n", + "plot(ids.edge_profiles.grid_ggd, n_e, colorbar_title=\"Electrons density / m^(-3)\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can overlap any grid on top of a quantity" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Choose backend\n", + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(ids.edge_profiles.grid_ggd, n_e) # Note default label in colorbar\n", + "plot!(space, GGDUtils.get_grid_subset_with_index(grid_ggd, 16), linecolor=:black, linewidth=2, linestyle=:solid, label=\"Separatix\", legend=true)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Adding equilibrium data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import SD4SOLPS: SD4SOLPS" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "eqdsk = \"samples/g002296.00200\"\n", + "SD4SOLPS.geqdsk_to_imas!(eqdsk, ids)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Extrapole core profile" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "SD4SOLPS.fill_in_extrapolated_core_profile!(ids, \"electrons.density\"; method=\"simple\")\n", + "SD4SOLPS.fill_in_extrapolated_core_profile!(ids, \"electrons.temperature\"; method=\"simple\")\n", + "# ... more profiles here as they become available in b2time" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Extend mesh outside SOLPS mesh to the device wall" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "SD4SOLPS.cached_mesh_extension!(ids, eqdsk, b2gmtry)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Loading a synthetic diagnostic" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using SynthDiag: SynthDiag" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Add interferometer chord details using a json file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "SynthDiag.add_interferometer!(\"samples/default_interferometer.json\", ids)" + ] + }, + { + "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 length and 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)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Load Langmuir Probes\n", + "\n", + "Same syntax as interferometer" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "SynthDiag.add_langmuir_probes!(\"samples/default_langmuir_probes.json\", ids)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " Data visualization recipes for langmuir probes have not been created yet and might not be created but you can still see them using built in plotting methodsan d knowledgeof IMAS data structure." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "gr() # Fast and can save pdf\n", + "# plotlyjs() # Use for interactive plot, can only save png\n", + "\n", + "plot(ids.langmuir_probes.embedded[1].time, ids.langmuir_probes.embedded[1].n_e.data, label=ids.langmuir_probes.embedded[1].name)\n", + "plot!(ylabel=\"Electron density / m^(-3)\", xlabel=\"Time / s\", legend=true)" + ] + } + ], + "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" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/example/makefile b/example/makefile new file mode 100644 index 0000000..8d6c7a2 --- /dev/null +++ b/example/makefile @@ -0,0 +1,25 @@ +SHELL := /bin/zsh +help: + @echo "Help Menu" + @echo + @echo "make env_with_cloned_repo: Creates a Julia environment with the cloned repositories" + @echo "make env_with_git_url: Creates a Julia environment with the git urls without creating local clones" + @echo "make clean: Deletes Project.toml and Manifest.toml for a fresh start" + @echo "make Clean: Deletes Project.toml, Manifest.toml, and any cloned repositories for a fresh start" + @echo + +env_with_cloned_repo: + chmod +x create_env_with_clone.sh + ./create_env_with_clone.sh + +env_with_git_url: + chmod +x create_env_with_git.sh + ./create_env_with_git.sh + +clean: + @echo "Deleting Project.toml and Manifest.toml" + rm Project.toml Manifest.toml + +Clean: clean + chmod +x rm_cloned_repo.sh + ./rm_cloned_repo.sh \ No newline at end of file diff --git a/example/rm_cloned_repo.sh b/example/rm_cloned_repo.sh new file mode 100644 index 0000000..351e0cc --- /dev/null +++ b/example/rm_cloned_repo.sh @@ -0,0 +1,32 @@ +#!/bin/bash + +echo "Deleting any cloned repo" +is_repo=$(git rev-parse --is-inside-work-tree 2>/dev/null) +if [[ ${is_repo} == "true" ]] +then + SD4SOLPS_path="$(git rev-parse --show-toplevel)" + project_path="$(dirname ${SD4SOLPS_path})" +else + SD4SOLPS_path="" + project_path="$(dirname $(pwd))" +fi + +echo "Removing extra cloned git repositories from "$project_path + +if [[ -z $SD4SOLPS_path ]] +then + # Check if $project_path"/SD4SOLPS.jl" exists + if [[ -d $project_path"/SD4SOLPS.jl" ]] + then + echo "Removing "$project_path"/SD4SOLPS.jl" + rm -rf $project_path"/SD4SOLPS.jl" + fi +fi +for repo in SynthDiag.jl SOLPS2IMAS.jl OMAS.jl GGDUtils.jl +do + if [[ -d $project_path"/"$repo ]] + then + echo "Removing "$project_path"/"$repo + rm -rf $project_path"/"$repo + fi +done \ No newline at end of file diff --git a/src/SD4SOLPS.jl b/src/SD4SOLPS.jl index 2fbeb29..0f87fe4 100644 --- a/src/SD4SOLPS.jl +++ b/src/SD4SOLPS.jl @@ -7,7 +7,7 @@ using Interpolations: Interpolations #import GGDUtils export find_files_in_allowed_folders -export geqdsk_to_imas +export geqdsk_to_imas! include("$(@__DIR__)/supersize_profile.jl") include("$(@__DIR__)/repair_eq.jl") @@ -63,12 +63,12 @@ function find_files_in_allowed_folders( end """ - geqdsk_to_imas() + geqdsk_to_imas!() Transfers the equilibrium reconstruction in an EFIT-style gEQDSK file into the IMAS DD structure. """ -function geqdsk_to_imas(eqdsk_file, dd; time_index=1) +function geqdsk_to_imas!(eqdsk_file, dd; time_index=1) # https://github.com/JuliaFusion/EFIT.jl/blob/master/src/io.jl g = EFIT.readg(eqdsk_file) # Copying ideas from OMFIT: omfit/omfit_classes/omfit_eqdsk.py / to_omas() @@ -122,6 +122,39 @@ function geqdsk_to_imas(eqdsk_file, dd; time_index=1) gq.q_min.rho_tor_norm = g.rhovn[qmin_idx] end + # X-points + xrs, xzs, xpsins, xseps = EFIT.x_points(g; within_limiter_only=false) + if length(xrs) > 0 + bx = eqt.boundary.x_point + resize!(bx, length(xrs)) + for i ∈ length(xrs) + bx[i].r = xrs[i] + bx[i].z = xzs[i] + end + nprim = sum(xseps .== 1) + if nprim > 0 + bsx = eqt.boundary_separatrix.x_point + resize!(bsx, nprim) + xrprim = xrs[xseps.==1] + xzprim = xzs[xseps.==1] + for i ∈ nprim + bsx[i].r = xrprim[i] + bsx[i].z = xzprim[i] + end + end + nsec = sum(xseps .== 2) + if nsec > 0 + bssx = eqt.boundary_secondary_separatrix.x_point + resize!(bssx, nsec) + xrsec = xrs[xseps.==2] + xzsec = xzs[xseps.==2] + for i ∈ nsec + bssx[i].r = xrsec[i] + bssx[i].z = xzsec[i] + end + end + end + # Boundary / LCFS eqt.boundary.outline.r = g.rbbbs eqt.boundary.outline.z = g.zbbbs @@ -134,7 +167,8 @@ function geqdsk_to_imas(eqdsk_file, dd; time_index=1) limiter.type.description = "first wall" resize!(limiter.unit, 1) limiter.unit[1].outline.r = g.rlim - return limiter.unit[1].outline.z = g.zlim + limiter.unit[1].outline.z = g.zlim + return end """ @@ -238,7 +272,7 @@ function preparation( println(" eqdsk = ", eqdsk) dd = SOLPS2IMAS.solps2imas(b2fgmtry, b2time, gridspec, b2mn) - geqdsk_to_imas(eqdsk, dd) + geqdsk_to_imas!(eqdsk, dd) println("Loaded input data into IMAS DD") fill_in_extrapolated_core_profile!(dd, "electrons.density"; method=core_method) @@ -246,6 +280,7 @@ function preparation( # ... more profiles here as they become available in b2time println("Extrapolated core profiles") + cached_mesh_extension!(dd, eqdsk_file, b2fgmtry) fill_in_extrapolated_edge_profile!(dd, "electrons.density"; method=core_method) # ... more profiles here println("Extrapolated edge profiles (but not really (placeholder only))") diff --git a/src/repair_eq.jl b/src/repair_eq.jl index d90b293..24b53bb 100644 --- a/src/repair_eq.jl +++ b/src/repair_eq.jl @@ -48,7 +48,7 @@ function add_rho_to_equilibrium!(dd::OMAS.dd) psi = eqt.profiles_1d.psi n = length(psi) if (length(eqt.profiles_1d.rho_tor_norm) > 0) - if max(eqt.profiles_1d.rho_tor_norm) > 0 + if maximum(eqt.profiles_1d.rho_tor_norm) > 0 println("Slice #", it, " already has a rho_tor_norm profile; skipping") continue end diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index 5e8cedc..cb75739 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -5,11 +5,16 @@ Utilities for extrapolating profiles # import CalculusWithJulia using OMAS: OMAS using Interpolations: Interpolations -using GGDUtils: GGDUtils +using GGDUtils: + GGDUtils, get_grid_subset_with_index, add_subset_element!, get_subset_boundary, + project_prop_on_subset!, get_subset_centers +using PolygonOps: PolygonOps +using JSON: JSON export extrapolate_core export fill_in_extrapolated_core_profile export mesh_psi_spacing +export find_x_points! """ cumul_integrate(x::AbstractVector, y::AbstractVector) @@ -132,9 +137,9 @@ function fill_in_extrapolated_core_profile!( grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] space = grid_ggd.space[space_idx] cell_subset = - SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 5) + get_grid_subset_with_index(grid_ggd, 5) midplane_subset = - SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 11) + get_grid_subset_with_index(grid_ggd, 11) if length(midplane_subset.element) < 1 throw( @@ -168,7 +173,7 @@ function fill_in_extrapolated_core_profile!( quantity_str = getproperty(quantity_str, Symbol(tag)) end - midplane_cell_centers, quantity = GGDUtils.project_prop_on_subset!( + midplane_cell_centers, quantity = project_prop_on_subset!( quantity_str, cell_subset, midplane_subset, @@ -205,10 +210,8 @@ function fill_in_extrapolated_core_profile!( in_bounds = (r .< maximum(r_eq)) .& (r .> minimum(r_eq)) .& (z .> minimum(z_eq)) .& (z .< maximum(z_eq)) - # println(in_bounds) psi_for_quantity = 10.0 .+ zeros(length(r)) psi_for_quantity[in_bounds] = rzpi.(r[in_bounds], z[in_bounds]) - # println(length(psi1_eq), ", ", length(rho1_eq)) rho_for_quantity = copy(psi_for_quantity) in_bounds = psi_for_quantity .<= 1.0 dpsi = diff(psi1_eq) @@ -279,6 +282,31 @@ function extrapolate_edge_exp( return y0 * exp(-x ./ lambda) end +""" + prep_flux_map() + +Reads equilibrium data and extracts/derives some useful quantities. +This is very basic, but it was being repeated and that's a no-no. +Returns: + + - R values of the equilibrium grid + - Z values of the eq grid + - normalized poloidal flux on the equilibrium grid + - a linear interpolation of norm pol flux vs. R and Z, ready to be evaluated +""" +function prep_flux_map(dd::OMAS.dd; eq_time_idx::Int64=1, eq_profiles_2d_idx::Int64=1) + eqt = dd.equilibrium.time_slice[eq_time_idx] + p2 = eqt.profiles_2d[eq_profiles_2d_idx] + r_eq = p2.grid.dim1 + z_eq = p2.grid.dim2 + psi = p2.psi + psia = eqt.global_quantities.psi_axis + psib = eqt.global_quantities.psi_boundary + psin_eq = (psi .- psia) ./ (psib - psia) + rzpi = Interpolations.LinearInterpolation((r_eq, z_eq), psin_eq) + return r_eq, z_eq, psin_eq, rzpi +end + #! format off """ mesh_psi_spacing( @@ -305,6 +333,10 @@ grid_ggd_idx: index of the grid_ggd to use. For a typical SOLPS run, the SOLPS g is used, then this index will need to be specified. space_idx: index of the space to use. For a typical SOLPS run, there will be only one space so this index will mostly remain at 1. +avoid_guard_cell: assume that the last cell is a guard cell so take end-2 and end-1 + instead of end and end-1 +spacing_rule: "edge" or "mean" to make spacing of new cells (in psi_N) be the same + as the spacing at the edge of the mesh, or the same as the average spacing """ #! format on function mesh_psi_spacing( @@ -313,6 +345,8 @@ function mesh_psi_spacing( eq_profiles_2d_idx::Int64=1, grid_ggd_idx::Int64=1, space_idx::Int64=1, + avoid_guard_cell::Bool=true, + spacing_rule="mean", ) # Inspect input if length(dd.equilibrium.time_slice) < eq_time_idx @@ -335,17 +369,7 @@ function mesh_psi_spacing( end # Get flux map - eqt = dd.equilibrium.time_slice[eq_time_idx] - p2 = eqt.profiles_2d[eq_profiles_2d_idx] - r_eq = p2.grid.dim1 - z_eq = p2.grid.dim2 - psi = p2.psi - psia = eqt.global_quantities.psi_axis - psib = eqt.global_quantities.psi_boundary - psin_eq = (psi .- psia) ./ (psib - psia) - rzpi = Interpolations.linear_interpolation((r_eq, z_eq), psin_eq) - println(minimum(r_eq), ", ", maximum(r_eq)) - println(minimum(z_eq), ", ", maximum(z_eq)) + r_eq, z_eq, psin_eq, rzpi = prep_flux_map(dd; eq_time_idx, eq_profiles_2d_idx) # Get a row of cells. Since the mesh should be aligned to the flux surfaces, # it shouldn't matter which row is used, although the divertor rows might be @@ -353,20 +377,612 @@ function mesh_psi_spacing( grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] space = grid_ggd.space[space_idx] midplane_subset = - SOLPS2IMAS.get_grid_subset_with_index(grid_ggd, 11) - midplane_cell_centers = GGDUtils.get_subset_centers(space, midplane_subset) + get_grid_subset_with_index(grid_ggd, 11) + midplane_cell_centers = get_subset_centers(space, midplane_subset) r_mesh = [midplane_cell_centers[i][1] for i ∈ eachindex(midplane_cell_centers)] z_mesh = [midplane_cell_centers[i][2] for i ∈ eachindex(midplane_cell_centers)] - println(minimum(r_mesh), ", ", maximum(r_mesh)) - println(minimum(z_mesh), ", ", maximum(z_mesh)) psin_mesh = rzpi.(r_mesh, z_mesh) # This should come out sorted, but with GGD, who knows. ii = sortperm(psin_mesh) psin = psin_mesh[ii] - dpsin = psin[end] - psin[end-1] + if spacing_rule == "edge" + if avoid_guard_cell + dpsin = psin[end-1] - psin[end-2] + else + dpsin = psin[end] - psin[end-1] + end + else + if avoid_guard_cell + dpsin = diff(psin[2:end-1]) + else + dpsin = diff(psin) + end + dpsin = sum(dpsin) / length(dpsin) + end return dpsin end +""" + pick_extension_psi_range() + +Defines the psi_N levels for an extended mesh. The range of psi_N levels starts +at the outer edge of the existing edge_profiles mesh at the midplane and goes +out to the most distant (in flux space) point on the limiting surface. +Returns a vector of psi_N levels. +""" +function pick_extension_psi_range( + dd::OMAS.dd; + eq_time_idx::Int64=1, + eq_profiles_2d_idx::Int64=1, + grid_ggd_idx::Int64=1, + space_idx::Int64=1, +) + r_eq, z_eq, psin_eq, rzpi = prep_flux_map(dd; eq_time_idx, eq_profiles_2d_idx) + + # Use wall to find maximum extent of contouring project + limiter = dd.wall.description_2d[1].limiter + wall_r = limiter.unit[1].outline.r + wall_z = limiter.unit[1].outline.z + wall_psin = rzpi.(wall_r, wall_z) + + # Use ggd mesh to find inner limit of contouring project + grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] + space = grid_ggd.space[space_idx] + midplane_subset = get_grid_subset_with_index(grid_ggd, 11) + midplane_cell_centers = get_subset_centers(space, midplane_subset) + psin_midplane = rzpi.(midplane_cell_centers[end][1], midplane_cell_centers[end][2]) + + # Choose contour levels + # The psi spacing function doesn't do much that's unique anymore. + # Should it be absorbed into this function intead? + dpsin = mesh_psi_spacing( + dd; + eq_time_idx=eq_time_idx, + eq_profiles_2d_idx=eq_profiles_2d_idx, + grid_ggd_idx=grid_ggd_idx, + space_idx=space_idx, + avoid_guard_cell=true, + spacing_rule="mean", + ) + lvlstart = maximum(psin_midplane * sign(dpsin)) / sign(dpsin) + dpsin + lvlend = maximum(wall_psin * sign(dpsin)) / sign(dpsin) + dpsin + nlvl = Int64(ceil((lvlend - lvlstart) / dpsin)) + psin_levels = collect(LinRange(lvlstart, lvlend, nlvl)) + # eqt = dd.equilibrium.time_slice[eq_time_idx] + # if hasproperty(eqt.boundary_secondary_separatrix, :psi) + # secondary_psi = eqt.boundary_secondary_separatrix.psi + # secondary_psin = (secondary_psi - psia) / (psib - psia) + # else + # secondary_psin = lvlend + dpsin + # end + return psin_levels +end + +""" + pick_mesh_ext_starting_points(dd; grid_ggd_idx, space_idx) + +Picks starting points for the radial lines of the mesh extension. The strategy +is to start from the outer edge of the existing mesh and follow the steepest +gradient (of psi_N) to extend these gridlines outward. +dd: a data dictionary instance with edge_profiles ggd and equilibrium loaded +grid_ggd_idx: index within ggd +space_idx: space number / index of the space to work with within edge_profiles +Returns a tuple with vectors of R and Z starting points. +""" +function pick_mesh_ext_starting_points(grid_ggd, space) + # Choose starting points for the orthogonal (to the contour) gridlines + # Use the existing cells of the standard mesh + all_cell_subset = get_grid_subset_with_index(grid_ggd, 5) + all_border_edges = get_subset_boundary(space, all_cell_subset) + core_edges = get_grid_subset_with_index(grid_ggd, 15) + outer_target = get_grid_subset_with_index(grid_ggd, 13) + inner_target = get_grid_subset_with_index(grid_ggd, 14) + ci = [ele.object[1].index for ele ∈ core_edges.element] + oi = [ele.object[1].index for ele ∈ outer_target.element] + ii = [ele.object[1].index for ele ∈ inner_target.element] + border_edges = [] + for i ∈ eachindex(all_border_edges) + bi = all_border_edges[i].object[1].index + if !(bi in oi) & !(bi in ii) & !(bi in ci) + border_edges = [border_edges; all_border_edges[i]] + end + end + + npol = length(border_edges) + r = zeros(npol) + z = zeros(npol) + corner_idx = 1 + for i ∈ 1:npol + edge_idx = border_edges[i].object[1].index + node_idx = space.objects_per_dimension[2].object[edge_idx].nodes[corner_idx] + geo = space.objects_per_dimension[1].object[node_idx].geometry + r[i, 1] = geo[1] + z[i, 1] = geo[2] + end + return r, z +end + +#!format off +""" + mesh_ext_follow_grad() + +Follows the steepest gradient from a set of starting points, dropping nodes at +approximately regular intervals in psi_N. Due to the numerical techniques used, the +node spacing may be imperfect (especially prone to error in regions where curvature +of psi_N is large compared to its gradient). +r_eq: Equilibrium reconstruction's grid, R coordinates +z_eq: Equilibrium reconstruction's grid, Z coordinates +psin_eq: Normalized poloidal flux in the equilibrium reconstruction as a function of R and Z +rstart: R coordinates of starting points for the gradient following. +zstart: Z coordinates of starting points +nlvl: number of nodes to drop while following the gradient +dpsin: node spacing in delta psi_N +rzpi: linear interpolation of psin_eq() as a function of r_eq and z_eq + This was probably already computed and I think time would be saved by reusing it. + If you don't already have it, you can pass in nothing and let this function calculate it. + +Returns two matrices with R and Z coordinates of the mesh extension +""" +#!format on +function mesh_ext_follow_grad( + r_eq::Vector{Float64}, + z_eq::Vector{Float64}, + psin_eq::Matrix, + rstart::Vector{Float64}, + zstart::Vector{Float64}, + nlvl::Int64, + dpsin::Float64, + rzpi=nothing, +) + npol = length(rstart) + mesh_r = zeros((npol, nlvl)) + mesh_z = zeros((npol, nlvl)) + for i ∈ 1:npol + mesh_r[i, 1] = rstart[i] + mesh_z[i, 1] = zstart[i] + end + + if rzpi === nothing + rzpi = Interpolations.LinearInterpolation((r_eq, z_eq), psin_eq) + end + + # Step along the paths of steepest descent to populate the mesh. + dpsindr, dpsindz = OMAS.gradient(r_eq, z_eq, psin_eq) + dpdr = Interpolations.linear_interpolation((r_eq, z_eq), dpsindr) + dpdz = Interpolations.linear_interpolation((r_eq, z_eq), dpsindz) + rlim = (minimum(r_eq), maximum(r_eq)) + zlim = (minimum(z_eq), maximum(z_eq)) + pfr = rzpi.(mesh_r[:, 1], mesh_z[:, 1]) .< 1 + for i ∈ 1:npol + if (mesh_r[i, 1] > rlim[1]) & (mesh_r[i, 1] < rlim[2]) & + (mesh_z[i, 1] > zlim[1]) & (mesh_z[i, 1] < zlim[2]) + # direction = rzpi(mesh_r[i, 1], mesh_z[i, 1]) >= 1 ? 1 : -1 + direction = pfr[i] ? -1 : 1 + for j ∈ 2:nlvl + # This is low resolution linear shooting that could go wrong + mesh_r[i, j] = mesh_r[i, j-1] + mesh_z[i, j] = mesh_z[i, j-1] + upscale = 5 # Should improve gradient following in regions of low gradient (near X-point) + for k ∈ 1:upscale + if (mesh_r[i, j] > rlim[1]) & (mesh_r[i, j] < rlim[2]) & + (mesh_z[i, j] > zlim[1]) & (mesh_z[i, j] < zlim[2]) + dpr = dpdr(mesh_r[i, j], mesh_z[i, j]) + dpz = dpdz(mesh_r[i, j], mesh_z[i, j]) + d = dpsin * direction / upscale / (dpr .^ 2 + dpz .^ 2) + mesh_r[i, j] += d * dpr + mesh_z[i, j] += d * dpz + end + end + end + end + # println("i=",i,"; r=",mesh_r[i, 1],":", mesh_r[i, end],",z=",mesh_z[i, 1], ":", mesh_z[i, end]) + end + return mesh_r, mesh_z +end + +""" + modify_mesh_ext_near_x! + +Modifies an extended mesh near a secondary X-point to compensate for the +tendency of the mesh to go nuts near the X-point. +eqt: equilibrium.time_slice information +mesh_r: matrix of R values for the extended mesh +mesh_z: matrix of Z values for the extended mesh +""" +function modify_mesh_ext_near_x!( + eqt::OMAS.equilibrium__time_slice, + mesh_r::Matrix{Float64}, + mesh_z::Matrix{Float64}, +) + # There's a special path; the one that probably should've gone through the + # secondary X-point (if there is one). Any numerical error will make this + # path miss the X-point and go off to somewhere crazy instead, so instead of + # that, let's draw a straight line from the edge of the SOLPS mesh to the + # X-point. + npol, nlvl = size(mesh_r) + bssx = eqt.boundary_secondary_separatrix.x_point + nsx = length(bssx) + println("there are ", nsx, " secondary x points") + if nsx > 0 + # There is hopefully only one X point on the secondary separatrix, but + # just in case some joker made a secondary snowflake or something crazy, + # we should pick which one we like best. + if nsx == 1 + xidx = 1 # Easy peasy + else + score = zeros(nsx) + rsx = [bssx[i].r for i ∈ 1:nsx] + zsx = [bssx[i].z for i ∈ 1:nsx] + # Being close to the vertically flipped position of the primary X-pt + # seems nice; let's reward that + bsx = eqt.boundary_separatrix.x_point + npx = length(bsx) + rpx = [bsx[i].r for i ∈ 1:npx] + zpx = [bsx[i].z for i ∈ 1:npx] + # Well, if we have a primary snowflake, there could be multiple + # primary X-points, so now we have to pick our favorite primary. + # Being close to the boundary is a good sign + min_ds2 = ones(npx) * 1e6 + for j ∈ 1:nsx + dr = [mesh_r[i, 1] - rpx[j] for i ∈ 1:npol] + dz = [mesh_z[i, 1] - zpx[j] for i ∈ 1:npol] + ds2 = dr .^ 2 .+ dz .^ 2 + min_ds2[j] = minimum(ds2) + end + primary_idx = argmin(min_ds2) + rpx = rpx[primary_idx] + zpx = zpx[primary_idx] + # There can be only one + + # Distance between secondary X-points and vert flip of favorite primary + ds2xx = (rsx .- rpx) .^ 2 .+ (zsx .- (-zpx)) .^ 2 + score += ds2xx * 5 + + # Being close to the boundary is good + min_ds2xb = ones(nsx) * 1e6 + for j ∈ 1:nsx + dr = [mesh_r[i, 1] - rsx[j] for i ∈ 1:npol] + dz = [mesh_z[i, 1] - zsx[j] for i ∈ 1:npol] + ds2 = dr .^ 2 .+ dz .^ 2 + min_ds2xb[j] = minimum(ds2) + end + score += min_ds2xb + # And the winner is the one with the lowest score + xidx = argmin(score) + end + # Pick which edges are closest to the X-point + rx = bssx[xidx].r + zx = bssx[xidx].z + ds2 = (mesh_r[:, 1] .- rx) .^ 2 .+ (mesh_z[:, 1] .- zx) .^ 2 + closest = argmin(ds2) + # Work out the mesh spacing + drm = mesh_r[closest, 2] - mesh_r[closest, 1] + dzm = mesh_z[closest, 2] - mesh_z[closest, 1] + dm = sqrt(drm^2 + dzm^2) + # Pick the angle + dr = rx - mesh_r[closest, 1] + dz = zx - mesh_z[closest, 1] + angle = atan(dz, dr) + # Replace the points + for i ∈ 2:nlvl + mesh_r[closest, i] = mesh_r[closest, 1] + dm * (i - 1) * cos(angle) + mesh_z[closest, i] = mesh_z[closest, 1] + dm * (i - 1) * sin(angle) + end + end +end + +#!format off +""" + record_regular_mesh!() + +Records arrays of mesh data from regular 2D arrays into the DD +grid_ggd: grid_ggd within edge_profiles +space: space in edge_profiles +mesh_r: Matrix of R values along a mesh. Should be 2D. The two dimensions are in + the radial and poloidal directions. +mesh_z: Z values to go with mesh_r. +cut: Poloidal index of a cut between two groups of disconnected cells. Poloidal + connections (faces, cells) will not be added between this poloidal index + and the next index. +""" +#!format on +function record_regular_mesh!( + grid_ggd, + space, + mesh_r::Matrix{Float64}, + mesh_z::Matrix{Float64}, + cut::Int64, +) + npol, nlvl = size(mesh_r) + o0 = space.objects_per_dimension[1] # Nodes + o1 = space.objects_per_dimension[2] # Edges + o2 = space.objects_per_dimension[3] # Cells (2D projection) + node_dim = 1 + edge_dim = 2 + cell_dim = 3 + + n_old_node = length(o0.object) + n_new_node = nlvl * npol + n_nodes = n_old_node + n_new_node + + n_old_edge = length(o1.object) + n_new_edge_i = nlvl * (npol - 2) + n_new_edge_j = (nlvl - 1) * npol + n_new_edge = n_new_edge_i + n_new_edge_j + n_edges = n_old_edge + n_new_edge + + n_old_cell = length(o2.object) + n_new_cell = (nlvl - 1) * (npol - 2) # -2 to account for disconnect @ pfr + n_cells = n_old_cell + n_new_cell + + # Define starting points and increments + node_start = n_old_node + 1 + edge_start1 = n_old_edge + 1 + edge_start2 = edge_start1 + n_new_edge_i + cell_start = n_old_cell + 1 + + n_per_i = nlvl + n_per_j = 1 + e1_per_i = nlvl # Edges along the npol direction / i direction + e2_per_i = nlvl - 1 # Edges along the nlvl direction / j direction + c_per_i = nlvl - 1 # # of cells < # of nodes + + # Get new subsets ready + n_existing_subsets = length(grid_ggd.grid_subset) + new_subset_indices = [-1, -2, -3, -4, -5, -201, -202, -203, -204, -205] + n_new_subsets = length(new_subset_indices) + resize!(grid_ggd.grid_subset, n_existing_subsets + n_new_subsets) + for i ∈ 1:n_new_subsets + grid_ggd.grid_subset[n_existing_subsets+i].identifier.index = + new_subset_indices[i] + end + ext_nodes_sub = get_grid_subset_with_index(grid_ggd, -201) + ext_edges_sub = get_grid_subset_with_index(grid_ggd, -202) + ext_xedges_sub = get_grid_subset_with_index(grid_ggd, -203) + ext_yedges_sub = get_grid_subset_with_index(grid_ggd, -204) + ext_cells_sub = get_grid_subset_with_index(grid_ggd, -205) + + # Preserve record of standard (non extended) mesh + for i ∈ 1:5 + std_sub = get_grid_subset_with_index(grid_ggd, -i) + orig_sub = get_grid_subset_with_index(grid_ggd, i) + resize!(std_sub.element, length(orig_sub.element)) + for j ∈ 1:length(orig_sub.element) + std_sub.element[j] = deepcopy(orig_sub.element[j]) + end + std_sub.identifier.index = -i + std_sub.dimension = deepcopy(orig_sub.dimension) + std_sub.metric = deepcopy(orig_sub.metric) + end + all_nodes_sub = get_grid_subset_with_index(grid_ggd, 1) + all_edges_sub = get_grid_subset_with_index(grid_ggd, 2) + all_xedges_sub = get_grid_subset_with_index(grid_ggd, 3) + all_yedges_sub = get_grid_subset_with_index(grid_ggd, 4) + all_cells_sub = get_grid_subset_with_index(grid_ggd, 5) + + nodes = resize!(o0.object, n_nodes) + edges = resize!(o1.object, n_edges) + cells = resize!(o2.object, n_cells) + space_idx = space.identifier.index + for i ∈ 1:npol + pastpc = i > cut + for j ∈ 1:nlvl + # Modified counters + ii = i - 1 # Offset due to indexing from 1 + jj = j - 1 # Offset due to indexing from 1 + iii = ii - 1 - pastpc # # of connections < than # of nodes, missing row @ PFR transition + jjj = jj - 1 # # of connections < # of nodes + + # Nodes + node_idx = node_start + ii * n_per_i + jj + nodes[node_idx].geometry = [mesh_r[i, j], mesh_z[i, j]] + add_subset_element!(ext_nodes_sub, space_idx, node_dim, node_idx) + add_subset_element!(all_nodes_sub, space_idx, node_dim, node_idx) + + # Edges + if (i > 1) & (i != cut) # i-1 to i in the npol direction + edge_idx1 = edge_start1 + iii * e1_per_i + jj + edges[edge_idx1].nodes = [node_idx, node_idx - n_per_i] + resize!(edges[edge_idx1].boundary, 2) + edges[edge_idx1].boundary[1].index = node_idx + edges[edge_idx1].boundary[2].index = node_idx - n_per_i + add_subset_element!(ext_edges_sub, space_idx, edge_dim, edge_idx1) + add_subset_element!(ext_xedges_sub, space_idx, edge_dim, edge_idx1) + add_subset_element!(all_edges_sub, space_idx, edge_dim, edge_idx1) + add_subset_element!(all_xedges_sub, space_idx, edge_dim, edge_idx1) + end + if (j > 1) # j-1 to j in the nlvl direction + edge_idx2 = edge_start2 + ii * e2_per_i + jjj + edges[edge_idx2].nodes = [node_idx, node_idx - n_per_j] + resize!(edges[edge_idx2].boundary, 2) + edges[edge_idx2].boundary[1].index = node_idx + edges[edge_idx2].boundary[2].index = node_idx - n_per_j + add_subset_element!(ext_edges_sub, space_idx, edge_dim, edge_idx2) + add_subset_element!(ext_yedges_sub, space_idx, edge_dim, edge_idx2) + add_subset_element!(all_edges_sub, space_idx, edge_dim, edge_idx2) + add_subset_element!(all_yedges_sub, space_idx, edge_dim, edge_idx2) + end + + # Cells + if (i > 1) & (i != cut) & (j > 1) + cell_idx = cell_start + iii * c_per_i + jjj + cells[cell_idx].nodes = [ + node_idx, + node_idx - n_per_i, + node_idx - n_per_i - n_per_j, + node_idx - n_per_j, + ] + resize!(cells[cell_idx].boundary, 4) + cells[cell_idx].boundary[1].index = edge_idx1 + cells[cell_idx].boundary[2].index = edge_idx2 + cells[cell_idx].boundary[3].index = edge_idx1 - 1 + cells[cell_idx].boundary[4].index = edge_idx2 - e2_per_i + + add_subset_element!(ext_cells_sub, space_idx, cell_dim, cell_idx) + add_subset_element!(all_cells_sub, space_idx, cell_dim, cell_idx) + end + end + end +end + +""" + convert_filename(filename::String) + +Converts a filename into a string that doesn't have illegal characters. +The main application is removing the path separator from source files with full +paths so the full paths* can be part of the new filename. This way, the input +files used to form some data can be part of the cache name, allowing quick lookup: +the cache filename is defined by the input files, and if it doesn't exist, it +needs to be generated. +""" +function convert_filename(filename::String) + filename_mod = replace(filename, "/" => "__") # Illegal on *nix bc it's the path separator + filename_mod = replace(filename_mod, ":" => "--") # Illegal on mac and windows + filename_mod = replace(filename_mod, "\\" => "__") # Illegal on windows + filename_mod = replace(filename_mod, "\"" => "''") # Illegal on windows + filename_mod = replace(filename_mod, "<" => "_lt_") # Illegal on windows + filename_mod = replace(filename_mod, ">" => "_gt_") # Illegal on windows + filename_mod = replace(filename_mod, "|" => "_pipe_") # Illegal on windows + filename_mod = replace(filename_mod, "?" => "_q_") # Illegal on windows + filename_mod = replace(filename_mod, "*" => "_a_") # Illegal on windows + return filename_mod +end + +#!format off +""" + cached_mesh_extension!() + +Adds an extended mesh to a data dictionary, possibly from a cached result. +dd: The data dictionary. It will be modified in place. +eqdsk_file: the name of the EQDSK file that was used to get equilibrium data in + the dd. +b2fgmtry: the name of the SOLPS geometry file that was used to get GGD info in + edge_profiles in the dd. +eq_time_idx: Index of the time slice in equilibrium +eq_profiles_2d_idx: Index of the 2D profile set in equilibrium + (there is usually only one) +grid_ggd_idx: Index of the grid_ggd set in edge_profiles +space_idx: Index of the space +clear_cache: delete any existing cache file (for use in testing) +""" +#!format on +function cached_mesh_extension!( + dd::OMAS.dd, + eqdsk_file::String, + b2fgmtry::String; + eq_time_idx::Int64=1, + eq_profiles_2d_idx::Int64=1, + grid_ggd_idx::Int64=1, + space_idx::Int64=1, + clear_cache=false, +) + path = "$(@__DIR__)/../data/" + eqdsk_file_mod = convert_filename(eqdsk_file) + b2fgmtry_mod = convert_filename(eqdsk_file) + cached_ext_name = path * eqdsk_file_mod * "_" * b2fgmtry_mod * ".mesh_ext.json" + if clear_cache + rm(cached_ext_name; force=true) + return cached_ext_name + end + if isfile(cached_ext_name) + # md = YAML.load_file(cached_ext_name) + md = JSON.parsefile(cached_ext_name) + pfr_transition = md["pfr_transition"] + mesh_r = convert(Matrix{Float64}, mapreduce(permutedims, vcat, md["r"])') + mesh_z = convert(Matrix{Float64}, mapreduce(permutedims, vcat, md["z"])') + npol = md["npol"] + nlvl = md["nlvl"] + record_regular_mesh!( + dd.edge_profiles.grid_ggd[grid_ggd_idx], + dd.edge_profiles.grid_ggd[grid_ggd_idx].space[space_idx], + mesh_r, + mesh_z, + pfr_transition, + ) + else + mesh_r, mesh_z, pfr_transition = mesh_extension_sol!( + dd; + eq_time_idx=eq_time_idx, + eq_profiles_2d_idx=eq_profiles_2d_idx, + grid_ggd_idx=grid_ggd_idx, + space_idx=space_idx, + ) + npol, nlvl = size(mesh_r) + data = Dict() + data["pfr_transition"] = pfr_transition + data["npol"] = npol + data["nlvl"] = nlvl + data["r"] = mesh_r + data["z"] = mesh_z + # YAML.write_file(cached_ext_name, data) + open(cached_ext_name, "w") do f + return JSON.print(f, data) + end + # fr = open("mesh_r.dat", "w") + # fz = open("mesh_z.dat", "w") + # for i ∈ 1:npol + # for j ∈ 1:nlvl + # print(fr, mesh_r[i, j], " ") + # print(fz, mesh_z[i, j], " ") + # end + # println(fr, "") + # println(fz, "") + # end + end + return cached_ext_name +end + +""" + function mesh_extension_sol() + +Extends the mesh out into the SOL +""" +function mesh_extension_sol!( + dd::OMAS.dd; + eq_time_idx::Int64=1, + eq_profiles_2d_idx::Int64=1, + grid_ggd_idx::Int64=1, + space_idx::Int64=1, +) + grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] + space = grid_ggd.space[space_idx] + eqt = dd.equilibrium.time_slice[eq_time_idx] + + r_eq, z_eq, psin_eq, rzpi = prep_flux_map(dd; eq_time_idx, eq_profiles_2d_idx) + psin_levels = pick_extension_psi_range( + dd; + eq_time_idx, + eq_profiles_2d_idx, + grid_ggd_idx, + space_idx, + ) + nlvl = length(psin_levels) + dpsin = psin_levels[2] - psin_levels[1] + grad_start_r, grad_start_z = pick_mesh_ext_starting_points(grid_ggd, space) + npol = length(grad_start_r) + mesh_r, mesh_z = mesh_ext_follow_grad( + r_eq, + z_eq, + psin_eq, + grad_start_r, + grad_start_z, + nlvl, + dpsin, + rzpi, + ) + modify_mesh_ext_near_x!(eqt, mesh_r, mesh_z) + + # Now we have all the nodes needed for the new mesh, but none are connected + # yet. They are, however, organized nicely in order, so it shouldn't be too + # hard. Any X-points in the domain should be aggressively ignored, so all the + # connections should be nice and simple. + # The PFR flag needs to be respected; pfr nodes don't connect to non-pfr nodes + pfr = rzpi.(mesh_r[:, 1], mesh_z[:, 1]) .< 1 + pfr_transition = argmax(abs.(diff(pfr))) + record_regular_mesh!(grid_ggd, space, mesh_r, mesh_z, pfr_transition) + return mesh_r, mesh_z, pfr_transition +end + """ fill_in_extrapolated_edge_profile!( dd::OMAS.dd, quantity_name::String; method::String="simple", diff --git a/test/runtests.jl b/test/runtests.jl index a8eebcf..ff6f825 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,6 +7,7 @@ using Test using Unitful: Unitful using Interpolations: Interpolations using ArgParse: ArgParse +using GGDUtils: GGDUtils, get_grid_subset_with_index function parse_commandline() s = ArgParse.ArgParseSettings(; description="Run tests. Default is all tests.") @@ -108,7 +109,7 @@ function define_default_sample_set() b2fgmtry, b2time, b2mn, gridspec, eqdsk = file_list eqdsk = splitdir(pathof(SD4SOLPS))[1] * - "/../sample/ITER_Lore_2296_00000/EQDSK/Baseline2008-li0.70.x4.mod2.eqdsk" + "/../sample/ITER_Lore_2296_00000/EQDSK/g002296.00200" return b2fgmtry, b2time, b2mn, gridspec, eqdsk end @@ -185,7 +186,7 @@ if args["core_profile_extension"] @test isfile(gridspec) @test isfile(eqdsk) dd = SOLPS2IMAS.solps2imas(b2fgmtry, b2time, gridspec, b2mn) - SD4SOLPS.geqdsk_to_imas(eqdsk, dd) + SD4SOLPS.geqdsk_to_imas!(eqdsk, dd) rho = dd.equilibrium.time_slice[1].profiles_1d.rho_tor_norm if !SD4SOLPS.check_rho_1d(dd; time_slice=1) @@ -221,10 +222,45 @@ if args["edge_profile_extension"] # Test for getting mesh spacing b2fgmtry, b2time, b2mn, gridspec, eqdsk = define_default_sample_set() dd = SOLPS2IMAS.solps2imas(b2fgmtry, b2time, gridspec, b2mn) - SD4SOLPS.geqdsk_to_imas(eqdsk, dd) + SD4SOLPS.geqdsk_to_imas!(eqdsk, dd) dpsin = SD4SOLPS.mesh_psi_spacing(dd) @test dpsin > 0.0 + # Extend the mesh + grid_ggd_idx = 1 + grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] + extended_subs = 1:5 + orig_subs = [ + deepcopy(get_grid_subset_with_index(grid_ggd, i)) for + i ∈ extended_subs + ] + cfn = SD4SOLPS.cached_mesh_extension!(dd, eqdsk, b2fgmtry; clear_cache=true) + println("cleared ext mesh cache: ", cfn) + SD4SOLPS.cached_mesh_extension!(dd, eqdsk, b2fgmtry; grid_ggd_idx=grid_ggd_idx) + for j ∈ extended_subs + orig_sub = orig_subs[j] + std_sub = get_grid_subset_with_index(grid_ggd, -j) + all_sub = get_grid_subset_with_index(grid_ggd, j) + ext_sub = get_grid_subset_with_index(grid_ggd, -200 - j) + orig_indices = [ele.object[1].index for ele ∈ orig_sub.element] + std_indices = [ele.object[1].index for ele ∈ std_sub.element] + all_indices = [ele.object[1].index for ele ∈ all_sub.element] + ext_indices = [ele.object[1].index for ele ∈ ext_sub.element] + @test std_sub.identifier.index == -j + @test all(orig_indices .== std_indices) + all_indices_reconstruct = [std_indices; ext_indices] + @test all(all_indices .== all_indices_reconstruct) + + # Verify that original and standard are separate and not refs to each other + arbitrary_change = 5 + arb_el = 2 + orig_sub.element[arb_el].object[1].index += arbitrary_change + @test orig_sub.element[arb_el].object[1].index != + std_sub.element[arb_el].object[1].index + orig_sub.element[arb_el].object[1].index -= arbitrary_change + end + + # Prepare profiles that need to be extended n_edge = 47 n_outer_prof = 13 quantity_edge = @@ -260,7 +296,7 @@ if args["heavy_utilities"] dd = OMAS.dd() eqdsk_file = splitdir(pathof(SD4SOLPS))[1] * "/../sample/geqdsk_iter_small_sample" - SD4SOLPS.geqdsk_to_imas(eqdsk_file, dd) + SD4SOLPS.geqdsk_to_imas!(eqdsk_file, dd) quantity = "electrons.density" prof_time_idx = eq_time_idx = 1 resize!(dd.core_profiles.profiles_1d, prof_time_idx) @@ -293,7 +329,7 @@ if args["repair_eq"] # Prepare sample dd = OMAS.dd() eqdsk = splitdir(pathof(SD4SOLPS))[1] * "/../sample/geqdsk_iter_small_sample" - SD4SOLPS.geqdsk_to_imas(eqdsk, dd) + SD4SOLPS.geqdsk_to_imas!(eqdsk, dd) # Make sure rho is missing nt = length(dd.equilibrium.time_slice) for it ∈ 1:nt @@ -324,7 +360,7 @@ if args["geqdsk_to_imas"] for sample_file ∈ sample_files println(sample_file) dd = OMAS.dd() - SD4SOLPS.geqdsk_to_imas(sample_file, dd; time_index=tslice) + SD4SOLPS.geqdsk_to_imas!(sample_file, dd; time_index=tslice) eqt = dd.equilibrium.time_slice[tslice] # global @@ -376,6 +412,16 @@ if args["geqdsk_to_imas"] # derived @test gq.q_axis == p1.q[1] + # X-points + bx = eqt.boundary.x_point + bsx = eqt.boundary_separatrix.x_point + bssx = eqt.boundary_secondary_separatrix.x_point + nxpt = length(bx) + nprim = length(bsx) + nsec = length(bssx) + @test nxpt >= 1 + @test nxpt >= (nprim + nsec) + # wall limiter = dd.wall.description_2d[1].limiter @test length(limiter.unit[1].outline.r) > 10