diff --git a/gpec_stability_scan.ipynb b/gpec_stability_scan.ipynb new file mode 100644 index 0000000..15b10bd --- /dev/null +++ b/gpec_stability_scan.ipynb @@ -0,0 +1,299 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "a1b2c3d4-0001-0001-0001-000000000001", + "metadata": {}, + "source": [ + "# Ideal MHD stability scan with ActorGPEC\n", + "\n", + "This notebook demonstrates how to use `FUSE.StudyDatabaseGenerator` together with `ActorGPEC` to scan a single equilibrium parameter and observe how the ideal-MHD stability energy δW(n=1) changes.\n", + "\n", + "**What we do:**\n", + "1. Load the D3D `:default` base case with boundary built from scalar parameters.\n", + "2. Vary the major radius R0 over ±12 % in 11 uniform steps using the `↔` notation.\n", + "3. Run `ActorGPEC` (n=1, no-wall vacuum) at every point using `StudyDatabaseGenerator` (serial execution).\n", + "4. Extract δW from `dd.mhd_linear` and plot it against R0.\n", + "\n", + "**Expected result:** a smooth curve — noisy/jumping values would signal a bug in the IMAS read path or the equilibrium solve." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a1b2c3d4-0001-0001-0001-000000000002", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "using FUSE\n", + "import IMAS\n", + "import JSON\n", + "using Plots\n", + "using Printf\n", + "FUSE.logging(Logging.Info; actors=Logging.Error);" + ] + }, + { + "cell_type": "markdown", + "id": "a1b2c3d4-0001-0001-0001-000000000003", + "metadata": {}, + "source": [ + "## 1. Base case and GPEC settings" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a1b2c3d4-0001-0001-0001-000000000004", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "ini, act = FUSE.case_parameters(:D3D, :default)\n", + "\n", + "# Build the boundary from scalar parameters (R0, a, κ, δ) so that\n", + "# varying ini.equilibrium.R0 actually changes the equilibrium geometry.\n", + "ini.equilibrium.boundary_from = :scalars\n", + "\n", + "# GPEC: n=1 toroidal mode, no conducting wall, vacuum perturbation only\n", + "act.ActorGPEC.nn_low = 1\n", + "act.ActorGPEC.nn_high = 1\n", + "act.ActorGPEC.vac_flag = true\n", + "act.ActorGPEC.wall_shape = \"nowall\"\n", + "act.ActorGPEC.verbose = false\n", + "act.ActorGPEC.psihigh = 0.95 # stay inside the FUSE equilibrium domain\n", + "act.ActorGPEC.qlow = 0.0 # no q threshold\n", + "\n", + "R0_base = ini.equilibrium.R0\n", + "println(\"Base major radius: R0 = $(round(R0_base; sigdigits=4)) m\")" + ] + }, + { + "cell_type": "markdown", + "id": "a1b2c3d4-0001-0001-0001-000000000005", + "metadata": {}, + "source": [ + "## 2. Parameter scan definition\n", + "\n", + "The `↔` operator marks a parameter as a scan variable. \n", + "`StudyDatabaseGenerator` samples `n_simulations` points uniformly between the two bounds." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a1b2c3d4-0001-0001-0001-000000000006", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Vary R0 ±12 % around the nominal value in 11 uniform steps\n", + "R0_min = 0.88 * R0_base\n", + "R0_max = 1.12 * R0_base\n", + "ini.equilibrium.R0 = R0_base ↔ [R0_min, R0_max]\n", + "\n", + "println(\"Scan range: R0 ∈ [$(round(R0_min; sigdigits=4)), $(round(R0_max; sigdigits=4))] m\")" + ] + }, + { + "cell_type": "markdown", + "id": "a1b2c3d4-0001-0001-0001-000000000007", + "metadata": {}, + "source": [ + "## 3. Study configuration and workflow" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a1b2c3d4-0001-0001-0001-000000000008", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "sty = FUSE.study_parameters(:DatabaseGenerator)\n", + "sty.n_workers = 0 # 0 = serial execution on the local machine\n", + "sty.n_simulations = 11\n", + "sty.save_folder = mktempdir() # temporary directory, change to a persistent path if desired\n", + "sty.file_save_mode = :overwrite\n", + "sty.database_policy = :separate_folders\n", + "\n", + "println(\"Results will be saved to: $(sty.save_folder)\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a1b2c3d4-0001-0001-0001-000000000009", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Workflow: initialise the equilibrium, then run ActorGPEC\n", + "function workflow_DatabaseGenerator(dd, ini, act)\n", + " FUSE.init(dd, ini, act)\n", + " FUSE.ActorGPEC(dd, act)\n", + " return nothing\n", + "end\n", + "\n", + "study = FUSE.StudyDatabaseGenerator(sty, ini, act)\n", + "study.workflow = workflow_DatabaseGenerator" + ] + }, + { + "cell_type": "markdown", + "id": "a1b2c3d4-0001-0001-0001-000000000010", + "metadata": {}, + "source": [ + "## 4. Run the scan" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a1b2c3d4-0001-0001-0001-000000000011", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "FUSE.run(study)" + ] + }, + { + "cell_type": "markdown", + "id": "a1b2c3d4-0001-0001-0001-000000000012", + "metadata": {}, + "source": [ + "## 5. Extract δW from saved results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a1b2c3d4-0001-0001-0001-000000000013", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "R0_values = Float64[]\n", + "dW_values = Float64[]\n", + "\n", + "for (k, folder) in enumerate(sort(filter(isdir, readdir(sty.save_folder; join=true))))\n", + " dd_file = joinpath(folder, \"dd.json\")\n", + " !isfile(dd_file) && (println(\" Case $k: dd.json not found\"); continue)\n", + " try\n", + " # Load only equilibrium + mhd_linear to avoid type-conversion issues\n", + " # from unrelated FUSE fields present in the full dd.json.\n", + " d = open(JSON.parse, dd_file)\n", + " case_dd = IMAS.dd()\n", + " IMAS.dict2imas(d[\"equilibrium\"], case_dd.equilibrium)\n", + " IMAS.dict2imas(d[\"mhd_linear\"], case_dd.mhd_linear)\n", + "\n", + " R0 = case_dd.equilibrium.time_slice[1].global_quantities.magnetic_axis.r\n", + "\n", + " if length(case_dd.mhd_linear.time_slice) > 0 &&\n", + " length(case_dd.mhd_linear.time_slice[1].toroidal_mode) > 0\n", + " δW = real(case_dd.mhd_linear.time_slice[1].toroidal_mode[1].energy_perturbed)\n", + " push!(R0_values, Float64(R0))\n", + " push!(dW_values, Float64(δW))\n", + " @printf \" Case %2d: R0 = %.3f m → δW(n=1) = %+.5f\\n\" k R0 δW\n", + " else\n", + " println(\" Case $k: no mhd_linear output (GPEC may have failed)\")\n", + " end\n", + " catch e\n", + " println(\" Case $k: failed to load — $(sprint(showerror, e))\")\n", + " end\n", + "end\n", + "\n", + "println(\"\\n$(length(dW_values))/$(sty.n_simulations) points succeeded.\")" + ] + }, + { + "cell_type": "markdown", + "id": "a1b2c3d4-0001-0001-0001-000000000014", + "metadata": {}, + "source": [ + "## 6. Plot δW vs R0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a1b2c3d4-0001-0001-0001-000000000015", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "length(dW_values) < 2 && error(\"Too few successful points to plot — check ActorGPEC setup.\")\n", + "\n", + "# Sort by R0 for a clean line plot\n", + "perm = sortperm(R0_values)\n", + "R0_sorted = R0_values[perm]\n", + "dW_sorted = dW_values[perm]\n", + "\n", + "p = plot(\n", + " R0_sorted, dW_sorted;\n", + " xlabel = \"Major radius R0 [m]\",\n", + " ylabel = \"δW(n=1) [normalized]\",\n", + " title = \"Ideal MHD stability vs major radius\\nD3D :default | n=1 | no-wall\",\n", + " marker = :circle,\n", + " markersize = 5,\n", + " lw = 2,\n", + " label = \"δW(n=1)\",\n", + " legend = :topright,\n", + ")\n", + "\n", + "hline!(p, [0.0];\n", + " linestyle = :dash,\n", + " color = :red,\n", + " lw = 1.5,\n", + " label = \"marginal stability (δW = 0)\",\n", + ")\n", + "\n", + "vline!(p, [R0_base];\n", + " linestyle = :dot,\n", + " color = :gray,\n", + " lw = 1.5,\n", + " label = \"nominal R0 = $(round(R0_base; sigdigits=4)) m\",\n", + ")\n", + "\n", + "display(p)" + ] + }, + { + "cell_type": "markdown", + "id": "a1b2c3d4-0001-0001-0001-000000000016", + "metadata": {}, + "source": [ + "**Reading the plot:**\n", + "- A smooth curve confirms that the IMAS read path and ActorGPEC are working correctly end-to-end.\n", + "- A noisy or jumping curve would indicate a problem in the spline construction inside `read_imas`.\n", + "- Where the curve crosses δW = 0 is the marginal-stability radius for the n=1 kink mode without a conducting wall." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia (8 threads) 1.11.1", + "language": "julia", + "name": "julia-_8-threads_-1.11" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.11.1" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}