Skip to content

Add Lao beta_p/li profile adjustment for VCs #79

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

Open
nicamo wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Add Lao beta_p/li profile adjustment for VCs #79

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
85de2f9
Add Lao beta_p/li profile adjustment for VCs
nicamo Jun 21, 2026
b5a6270
Correct profile-adjusted virtual circuits
nicamo Aug 4, 2026
e1607f2
Merge remote-tracking branch 'origin/main' into lao-betap-li-static-p…
nicamo Aug 4, 2026
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
249 changes: 249 additions & 0 deletions examples/example09 - virtual_circuits_MASTU.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -575,6 +575,255 @@
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Virtual circuits at fixed beta_p and li\n",
"\n",
"The VCs above are calculated at fixed profile parameters: when each coil current is perturbed, the same profile object is used again in the static GS solve. In some workflows it is more natural to hold global profile quantities fixed instead. Here we show how to build VCs at fixed `beta_p` and `li`.\n",
"\n",
"This route uses a `Lao85` profile with two alpha coefficients and two beta coefficients. The total plasma current is still enforced by the usual Lao global scaling, while the Lao coefficients are refitted after each coil-current perturbation so that:\n",
"\n",
"- `eq.poloidalBeta1()` stays fixed. This is the same beta_p definition used by `ConstrainBetapIp`.\n",
"- `eq.internalInductance2()` stays fixed. This is the EFIT/A-EQDSK-style `li` convention used here.\n",
"\n",
"The profile refit is local, so it is best suited to VC perturbations close to the reference equilibrium. By default the first refit builds the local Jacobian of `[beta_p, li]` with respect to `[alpha0, alpha1, beta0, beta1]`, and subsequent refits reuse that Jacobian for speed."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-21T16:03:58.789930Z",
"iopub.status.busy": "2026-06-21T16:03:58.789785Z",
"iopub.status.idle": "2026-06-21T16:03:59.998170Z",
"shell.execute_reply": "2026-06-21T16:03:59.997853Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Fixed beta_p target = 0.185358\n",
"Fixed li target = 0.910044\n",
"Initial Lao alpha coefficients = [ 1.04290355 -0.50239287]\n",
"Initial Lao beta coefficients = [ 2.85106681e-06 -1.37343059e-06]\n"
]
}
],
"source": [
"from freegsnke.jtor_update import Lao85\n",
"from freegsnke.virtual_circuits import make_lao85_betap_li_profile_adjuster\n",
"\n",
"# Build a two-alpha/two-beta Lao85 profile close to the ConstrainPaxisIp profile\n",
"# used above. The optional Lao logic terms are appended by Lao85 itself.\n",
"alpha_lao, beta_lao = profiles.Lao_parameters(\n",
" n_alpha=2,\n",
" n_beta=2,\n",
" alpha_logic=True,\n",
" beta_logic=True,\n",
")\n",
"\n",
"# Work on a copy so the rest of the notebook can continue to use `eq` and `profiles`.\n",
"eq_lao = eq.create_auxiliary_equilibrium()\n",
"profiles_lao = Lao85(\n",
" eq=eq_lao,\n",
" Ip=profiles.Ip,\n",
" fvac=profiles.fvac(),\n",
" alpha=alpha_lao,\n",
" beta=beta_lao,\n",
" alpha_logic=True,\n",
" beta_logic=True,\n",
" Ip_logic=True,\n",
")\n",
"\n",
"GSStaticSolver.forward_solve(\n",
" eq_lao,\n",
" profiles_lao,\n",
" target_relative_tolerance=1e-5,\n",
" suppress=True,\n",
")\n",
"\n",
"fixed_betap = eq_lao.poloidalBeta1()\n",
"fixed_li = eq_lao.internalInductance2()\n",
"print(f\"Fixed beta_p target = {fixed_betap:.6g}\")\n",
"print(f\"Fixed li target = {fixed_li:.6g}\")\n",
"print(f\"Initial Lao alpha coefficients = {alpha_lao}\")\n",
"print(f\"Initial Lao beta coefficients = {beta_lao}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `profile_adjuster` below is the object that changes the VC behaviour. It is called before the baseline solve and after each coil-current perturbation, before the target values are evaluated. In this example it refits the `Lao85` coefficients to keep `beta_p` and `li` fixed.\n",
"\n",
"The `x_scale` entry is useful because the pressure and FF' Lao coefficients can have very different magnitudes. The reusable metric Jacobian is stored on the adjuster after the first refit."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-21T16:03:59.999663Z",
"iopub.status.busy": "2026-06-21T16:03:59.999550Z",
"iopub.status.idle": "2026-06-21T16:04:00.001962Z",
"shell.execute_reply": "2026-06-21T16:04:00.001615Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Metric Jacobian cached before first refit: False\n"
]
}
],
"source": [
"lao_coefficient_scale = np.maximum(np.abs(np.r_[alpha_lao, beta_lao]), 1.0)\n",
"\n",
"fixed_betap_li_adjuster = make_lao85_betap_li_profile_adjuster(\n",
" reference_eq=eq_lao,\n",
" betap=fixed_betap,\n",
" li=fixed_li,\n",
" li_method=\"internalInductance2\",\n",
" optimizer_kwargs={\n",
" \"x_scale\": lao_coefficient_scale,\n",
" \"max_nfev\": 12,\n",
" },\n",
" use_metric_jacobian=True,\n",
" reuse_metric_jacobian=True,\n",
")\n",
"\n",
"print(\"Metric Jacobian cached before first refit:\", fixed_betap_li_adjuster.reusable_metric_jacobian[0] is not None)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we calculate a VC in the same way as before, but pass the adjuster through `profile_adjuster`. The resulting shape matrix is therefore calculated at fixed `beta_p` and fixed `li`, rather than at fixed Lao coefficients.\n",
"\n",
"This short worked example uses the robust midplane targets `Rin` and `Rout`. The same `profile_adjuster` route can be used with other target calculators, but x-point targets require each finite-difference perturbation to preserve an identifiable x-point. Each coil-current perturbation starts from the same adjusted baseline profile, and the baseline targets are evaluated only after that profile has been fitted."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def plasma_midplane_descriptors(eq):\n",
" return np.array(eq.innerOuterSeparatrix())\n",
"\n",
"\n",
"coils_fixed_betap_li = ['PX', 'D1', 'D2', 'D3']\n",
"\n",
"target_names_fixed_betap_li = [\"Rin\", \"Rout\"]\n",
"\n",
"VCs.calculate_VC(\n",
" eq=eq_lao,\n",
" profiles=profiles_lao,\n",
" coils=coils_fixed_betap_li,\n",
" target_names=target_names_fixed_betap_li,\n",
" target_calculator=plasma_midplane_descriptors,\n",
" starting_dI=None,\n",
" min_starting_dI=50,\n",
" verbose=True,\n",
" name=\"VC_for_lower_targets_fixed_betap_li\",\n",
" profile_adjuster=fixed_betap_li_adjuster,\n",
")\n",
"\n",
"print(\"Metric Jacobian cached:\", fixed_betap_li_adjuster.reusable_metric_jacobian[0] is not None)\n",
"print(\"Number of profile refits:\", len(fixed_betap_li_adjuster.fit_results))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Applying the VC works as before. Since the VC object stores the profile adjuster, `apply_VC` will also refit the Lao coefficients after applying the requested coil-current shifts."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"requested_target_shifts_fixed_betap_li = [0.02, -0.02]\n",
"\n",
"eq_lao_new, profiles_lao_new, new_target_values_fixed_betap_li, old_target_values_fixed_betap_li = VCs.apply_VC(\n",
" eq=eq_lao,\n",
" profiles=profiles_lao,\n",
" VC_object=VCs.VC_for_lower_targets_fixed_betap_li,\n",
" requested_target_shifts=requested_target_shifts_fixed_betap_li,\n",
" verbose=True,\n",
")\n",
"\n",
"print(\"Requested target shifts:\", requested_target_shifts_fixed_betap_li)\n",
"print(\"Actual target shifts: \", new_target_values_fixed_betap_li - old_target_values_fixed_betap_li)\n",
"print()\n",
"print(f\"Original beta_p = {eq_lao.poloidalBeta1():.6g}\")\n",
"print(f\"New beta_p = {eq_lao_new.poloidalBeta1():.6g}\")\n",
"print(f\"Original li = {eq_lao.internalInductance2():.6g}\")\n",
"print(f\"New li = {eq_lao_new.internalInductance2():.6g}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Visualise that the requested VC move is applied while beta_p and li remain fixed.\n",
"metric_names = [r\"$\\beta_p$\", r\"$l_i$\"]\n",
"old_metrics = np.array([eq_lao.poloidalBeta1(), eq_lao.internalInductance2()])\n",
"new_metrics = np.array([eq_lao_new.poloidalBeta1(), eq_lao_new.internalInductance2()])\n",
"\n",
"fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 9), dpi=80)\n",
"\n",
"ax1.grid(True, which=\"both\", alpha=0.5)\n",
"ax1.scatter(\n",
" target_names_fixed_betap_li,\n",
" requested_target_shifts_fixed_betap_li,\n",
" color=\"red\",\n",
" marker=\"o\",\n",
" s=150,\n",
" label=\"Requested\",\n",
")\n",
"ax1.scatter(\n",
" target_names_fixed_betap_li,\n",
" new_target_values_fixed_betap_li - old_target_values_fixed_betap_li,\n",
" color=\"royalblue\",\n",
" marker=\"o\",\n",
" s=75,\n",
" label=\"Actual\",\n",
")\n",
"ax1.set_xlabel(\"Target\")\n",
"ax1.set_ylabel(\"Shift [m]\")\n",
"ax1.legend()\n",
"\n",
"x = np.arange(len(metric_names))\n",
"ax2.grid(True, which=\"both\", alpha=0.5)\n",
"ax2.bar(x - 0.18, old_metrics, width=0.36, label=\"Before\")\n",
"ax2.bar(x + 0.18, new_metrics, width=0.36, label=\"After\")\n",
"ax2.set_xticks(x, metric_names)\n",
"ax2.set_ylabel(\"Value\")\n",
"ax2.legend()\n",
"\n",
"for i, (old_value, new_value) in enumerate(zip(old_metrics, new_metrics)):\n",
" rel_change = (new_value - old_value) / old_value\n",
" ax2.text(i, max(old_value, new_value), f\"{rel_change:+.2e}\", ha=\"center\", va=\"bottom\")\n",
"\n",
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down
3 changes: 0 additions & 3 deletions freegsnke/equilibrium_update.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -279,9 +279,6 @@ def psi_func(self, R, Z, *args, **kwargs):
> 1e-5
)
if check:
print(
"Dicrepancy between psi_func and plasma_psi detected. psi_func has been re-set."
)
# redefine interpolating function
self.psi_func_interp = interpolate.RectBivariateSpline(
self.R[:, 0], self.Z[0, :], self.plasma_psi
Expand Down
Loading
Loading