Skip to content

Commit b9ad1f2

Browse files
kpentlandkpentland
authored
Merge pull request #60 from timothy-nunn/weight-inverse-solver
Inverse Solver: Allow different types of constraints to be weighted in the inverse solver (can also now weight different isoflux constraints individually).
2 parents 35db8b8 + e430c0c commit b9ad1f2

2 files changed

Lines changed: 294 additions & 42 deletions

File tree

‎examples/example01b - advanced_static_inverse_solve.ipynb‎

Lines changed: 225 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -55,11 +55,22 @@
5555
"GSStaticSolver = GSstaticsolver.NKGSsolver(eq) "
5656
]
5757
},
58+
{
59+
"cell_type": "code",
60+
"execution_count": null,
61+
"metadata": {},
62+
"outputs": [],
63+
"source": [
64+
"# as before, let's fix the Solenoid current\n",
65+
"eq.tokamak.set_coil_current('Solenoid', 5000)\n",
66+
"eq.tokamak['Solenoid'].control = False # ensures the current in the Solenoid is fixed"
67+
]
68+
},
5869
{
5970
"cell_type": "markdown",
6071
"metadata": {},
6172
"source": [
62-
"### Constraints\n",
73+
"### (Normalised) poloidal flux constraints at specific locations\n",
6374
"\n",
6475
"In addition to the `null points`, `isoflux set`, and `coil_current_limits` constraints introduced in the previous notebook, additional methods are available for more advanced control of the inverse problem.\n",
6576
"\n",
@@ -81,6 +92,10 @@
8192
" $$\n",
8293
" \\hat{\\psi}(R_j, Z_j) \\leq \\hat{\\psi}_j^{\\text{target}}.\n",
8394
" $$\n",
95+
" or\n",
96+
" $$\n",
97+
" \\hat{\\psi}(R_j, Z_j) \\geq \\hat{\\psi}_j^{\\text{target}}.\n",
98+
" $$\n",
8499
" Note that $\\psi_{axis}$ and $\\psi_{boundary}$ are the values of the flux on the magnetic axis and plasma boundary, respectively.\n",
85100
" As we'll see, this can be useful for setting explicit constraints on flux behaviour near the wall if there are certain safety limits to adhere to.\n",
86101
"\n",
@@ -91,17 +106,6 @@
91106
"In the following, we will show how to use the normalised flux constraints to control the normalised flux on the divertor nose."
92107
]
93108
},
94-
{
95-
"cell_type": "code",
96-
"execution_count": null,
97-
"metadata": {},
98-
"outputs": [],
99-
"source": [
100-
"# as before, let's fix the Solenoid current\n",
101-
"eq.tokamak.set_coil_current('Solenoid', 5000)\n",
102-
"eq.tokamak['Solenoid'].control = False # ensures the current in the Solenoid is fixed"
103-
]
104-
},
105109
{
106110
"cell_type": "code",
107111
"execution_count": null,
@@ -119,12 +123,15 @@
119123
"# set desired null_points locations (this can include X-point and O-point locations)\n",
120124
"null_points = [[Rx, Rx], [Zx, -Zx]]\n",
121125
"\n",
122-
"# set desired isoflux constraints with format \n",
126+
"# Set desired isoflux constraints with format \n",
123127
"# isoflux_set = [isoflux_0, isoflux_1 ... ] \n",
124-
"# with each isoflux_i = [R_coords, Z_coords]\n",
128+
"# with each isoflux_i = [R_coords, Z_coords, weights]\n",
125129
"isoflux_set = np.array([\n",
126-
" [[Rx, Rx, Rin, Rout], [Zx, -Zx, 0.,0.]]\n",
127-
" ])\n",
130+
" [\n",
131+
" [Rx, Rx, Rin, Rout], \n",
132+
" [Zx, -Zx, 0., 0.],\n",
133+
" ]\n",
134+
"])\n",
128135
"\n",
129136
"# set the coil current limits (upper and lower)\n",
130137
"# coil ordering in this case: PX, D1, D2, D3, Dp, D5, D6, D7, P4, P5, P6\n",
@@ -147,16 +154,9 @@
147154
" isoflux_set=isoflux_set,\n",
148155
" coil_current_limits=coil_current_limits,\n",
149156
" psi_norm_limits=psi_norm_limits,\n",
157+
" mu_psi_norm = 1e7, # penalise how much the normalised psi constraint is violated (default 1e6)\n",
150158
")\n",
151-
"\n",
152-
"# if you find coil limits are being violated or an undesireable solution is being produced,\n",
153-
"# you can try increasing the penalty factor for violating the coil limits\n",
154-
"# (here we just set it to its default value of 1e5)\n",
155-
"constrain.mu_coils = 1e5\n",
156-
"\n",
157-
"# Similarly, you can define how much normalised psi constraints are penalised\n",
158-
"# (the default value is 1e6)\n",
159-
"constrain.mu_psi_norm = 1e7"
159+
"\n"
160160
]
161161
},
162162
{
@@ -205,7 +205,7 @@
205205
"cell_type": "markdown",
206206
"metadata": {},
207207
"source": [
208-
"Finally, you can check that the normalised $\\psi$ constraint has indeed been satisfied."
208+
"We can now check that the normalised $\\psi$ constraint has indeed been satisfied!"
209209
]
210210
},
211211
{
@@ -219,6 +219,205 @@
219219
" print(f\"Psi norm at ({psi_con[0]}, {psi_con[1]}) = {eq.psiNRZ(psi_con[0], psi_con[1]):.3f} ({sign} {psi_con[2]})\")"
220220
]
221221
},
222+
{
223+
"cell_type": "markdown",
224+
"metadata": {},
225+
"source": [
226+
"#### Weighting constraints within the solver\n",
227+
"\n",
228+
"It is possible to specify **weights** on the different types of constraints within the inverse solver, enabling the solver to prioritise certain constraints over others.\n",
229+
"\n",
230+
"There are two ways of doing this:\n",
231+
"\n",
232+
"1. **Weighting different types of constraints relative to one another:**\n",
233+
"\n",
234+
" Different classes of constraint can be weighted relative to each other (e.g. null points, isoflux, and psi constraints). This allows the solver to be instructed that one type of constraint is more important than others. For example, a user may want null point constraints to be strictly satisfied, but may be less concerned about the isoflux set being met with the same strictness.\n",
235+
"\n",
236+
"2. **Weighting individual isoflux constraints within a set:**\n",
237+
"\n",
238+
" Individual isoflux constraints within a set can optionally be weighted relative to one another, instructing the solver to prioritise satisfying some subset over others.\n",
239+
" \n",
240+
" For example, consider a set of three isoflux constraints: one at the X-point, one at the outer midplane, and one at a desired strike point. The strike point constraint could be assigned a lower weight to allow the solver to find a solution with a strike point close to — but not necessarily exactly on — the target location."
241+
]
242+
},
243+
{
244+
"cell_type": "markdown",
245+
"metadata": {},
246+
"source": [
247+
"#### Option 1: weighting different types of constraints\n",
248+
"\n",
249+
"Here, we'll show how to weight the null point constraints **more** than the isoflux constraints. This solver then prioritises fitting the null points exactly and the isoflux points less so. "
250+
]
251+
},
252+
{
253+
"cell_type": "code",
254+
"execution_count": null,
255+
"metadata": {},
256+
"outputs": [],
257+
"source": [
258+
"# core constraints\n",
259+
"Rx = 0.55 # X-point radius\n",
260+
"Zx = 1.2 # X-point height\n",
261+
"Rout = 1.4 # outboard midplane radius\n",
262+
"Rin = 0.34 # inboard midplane radius\n",
263+
"\n",
264+
"# set desired null_points locations (this can include X-point and O-point locations)\n",
265+
"null_points = [[Rx, Rx], [Zx, -Zx]]\n",
266+
"\n",
267+
"# set desired isoflux constraints with format \n",
268+
"# isoflux_set = [isoflux_0, isoflux_1 ... ] \n",
269+
"# with each isoflux_i = [R_coords, Z_coords]\n",
270+
"isoflux_set = np.array([\n",
271+
" [\n",
272+
" [Rx, Rx, Rin, Rout, 0.75, 1.0, 0.75, 1.0], \n",
273+
" [Zx, -Zx, 0.0, 0.0, -1.6, -2.0, 1.6, 2.0],\n",
274+
" ]\n",
275+
"])\n",
276+
"\n",
277+
"# instantiate the freegsnke constrain object\n",
278+
"constrain = Inverse_optimizer(\n",
279+
" null_points=null_points,\n",
280+
" isoflux_set=isoflux_set,\n",
281+
" coil_current_limits=coil_current_limits,\n",
282+
" weight_isoflux=0.2, # <-- weight the isofluxes less than the null points (default 1.0)\n",
283+
" weight_nulls=1.0, # <-- weight the null points more than the null points (default 1.0)\n",
284+
" # weight_psi=1.0 # <-- weight for psi values not used here\n",
285+
")\n",
286+
"\n"
287+
]
288+
},
289+
{
290+
"cell_type": "code",
291+
"execution_count": null,
292+
"metadata": {},
293+
"outputs": [],
294+
"source": [
295+
"# solve!\n",
296+
"GSStaticSolver.solve(eq=eq, \n",
297+
" profiles=profiles, \n",
298+
" constrain=constrain, \n",
299+
" target_relative_tolerance=1e-6,\n",
300+
" target_relative_psit_update=1e-3,\n",
301+
" verbose=True,\n",
302+
" l2_reg=np.array([1e-12]*10+[1e-6]), \n",
303+
" )"
304+
]
305+
},
306+
{
307+
"cell_type": "markdown",
308+
"metadata": {},
309+
"source": [
310+
"Notice how the null points constraints are satisfied pefectly while the isoflux constraints are \"less so\" (especially in the divertor region)."
311+
]
312+
},
313+
{
314+
"cell_type": "code",
315+
"execution_count": null,
316+
"metadata": {},
317+
"outputs": [],
318+
"source": [
319+
"# plot the resulting equilbria \n",
320+
"fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n",
321+
"ax1.grid(True, which='both')\n",
322+
"eq.plot(axis=ax1, show=False)\n",
323+
"eq.tokamak.plot(axis=ax1, show=False)\n",
324+
"constrain.plot(axis=ax1,show=True)\n",
325+
"ax1.set_xlim(0.1, 2.15)\n",
326+
"ax1.set_ylim(-2.25, 2.25)\n",
327+
"plt.tight_layout()"
328+
]
329+
},
330+
{
331+
"cell_type": "markdown",
332+
"metadata": {},
333+
"source": [
334+
"#### Option 2: weighting different isoflux constraints\n",
335+
"\n",
336+
"Here, we'll show how to weight individual isoflux constraints differently from one another. For example, here we want to prioritise the core shape over the exact location of the strikepoint. "
337+
]
338+
},
339+
{
340+
"cell_type": "code",
341+
"execution_count": null,
342+
"metadata": {},
343+
"outputs": [],
344+
"source": [
345+
"# core constraints\n",
346+
"Rx = 0.55 # X-point radius\n",
347+
"Zx = 1.2 # X-point height\n",
348+
"Rout = 1.4 # outboard midplane radius\n",
349+
"Rin = 0.34 # inboard midplane radius\n",
350+
"\n",
351+
"# set desired null_points locations (this can include X-point and O-point locations)\n",
352+
"null_points = [[Rx, Rx], [Zx, -Zx]]\n",
353+
"\n",
354+
"# set desired isoflux constraints with format \n",
355+
"# isoflux_set = [isoflux_0, isoflux_1 ... ] \n",
356+
"# with each isoflux_i = [R_coords, Z_coords, weights]\n",
357+
"isoflux_set = np.array([\n",
358+
" [\n",
359+
" [Rx, Rx, Rin, Rout, 1.0, 1.0, 0.75, 1.0, 0.75, 1.0], \n",
360+
" [Zx, -Zx, 0.0, 0.0, -0.9, 0.9, -1.6, -2.0, 1.6, 2.0],\n",
361+
" [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.1, 1.0, 0.1],\n",
362+
" ]\n",
363+
"])\n",
364+
"\n",
365+
"# instantiate the freegsnke constrain object\n",
366+
"constrain = Inverse_optimizer(\n",
367+
" null_points=null_points,\n",
368+
" isoflux_set=isoflux_set,\n",
369+
" coil_current_limits=coil_current_limits,\n",
370+
")"
371+
]
372+
},
373+
{
374+
"cell_type": "code",
375+
"execution_count": null,
376+
"metadata": {},
377+
"outputs": [],
378+
"source": [
379+
"# solve!\n",
380+
"GSStaticSolver.solve(eq=eq, \n",
381+
" profiles=profiles, \n",
382+
" constrain=constrain, \n",
383+
" target_relative_tolerance=1e-6,\n",
384+
" target_relative_psit_update=1e-3,\n",
385+
" verbose=True, # print output\n",
386+
" l2_reg=np.array([1e-12]*10+[1e-6]), \n",
387+
" )"
388+
]
389+
},
390+
{
391+
"cell_type": "markdown",
392+
"metadata": {},
393+
"source": [
394+
"Notice how the strikepoint isoflux constraint is no longer satisfied."
395+
]
396+
},
397+
{
398+
"cell_type": "code",
399+
"execution_count": null,
400+
"metadata": {},
401+
"outputs": [],
402+
"source": [
403+
"# plot the resulting equilbria \n",
404+
"fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n",
405+
"ax1.grid(True, which='both')\n",
406+
"eq.plot(axis=ax1, show=False)\n",
407+
"eq.tokamak.plot(axis=ax1, show=False)\n",
408+
"constrain.plot(axis=ax1,show=True)\n",
409+
"ax1.set_xlim(0.1, 2.15)\n",
410+
"ax1.set_ylim(-2.25, 2.25)\n",
411+
"plt.tight_layout()"
412+
]
413+
},
414+
{
415+
"cell_type": "code",
416+
"execution_count": null,
417+
"metadata": {},
418+
"outputs": [],
419+
"source": []
420+
},
222421
{
223422
"cell_type": "code",
224423
"execution_count": null,

0 commit comments

Comments
 (0)