Skip to content

Add ActorGPEC ideal MHD stability scan notebook #18

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
calin1989 wants to merge 1 commit into
base: master
Choose a base branch
from
Draft

Add ActorGPEC ideal MHD stability scan notebook #18

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
299 changes: 299 additions & 0 deletions gpec_stability_scan.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
}