Inverse solver: enabling second order magnetic null constraints (for snowflake divertors)#56
Open
kpentland wants to merge 13 commits into
Open
Inverse solver: enabling second order magnetic null constraints (for snowflake divertors)#56kpentland wants to merge 13 commits into
kpentland wants to merge 13 commits into
Conversation
kpentland
added
ready-for-final-tests
kpentland
added
ready-for-final-tests
Contributor
|
Sorry @kpentland I forgot about this, I have put it on my list for review tomorrow! |
timothy-nunn
suggested changes
Apr 15, 2026
timothy-nunn
suggested changes
May 6, 2026
Comment on lines
+104
to
+110
| isoflux_set : list or np.array, optional | ||
| list of isoflux objects, each with structure | ||
| [Rcoords, Zcoords] | ||
| with Rcoords and Zcoords being 1D lists of the coords of all points that are requested to have the same flux value | ||
| null_points : list or np.array, optional | ||
| structure [Rcoords, Zcoords], with Rcoords and Zcoords being 1D lists | ||
| Sets the coordinates of the desired null points, including both Xpoints and Opoints |
Contributor
There was a problem hiding this comment.
isoflux_set and null_points is already documented above, is this a duplicate?
Comment on lines
+115
to
+118
| psi_vals : list or np.array, optional | ||
| structure [Rcoords, Zcoords, psi_values] | ||
| with Rcoords, Zcoords and psi_values having the same shape | ||
| Sets the desired values of psi for a set of coordinates, possibly an entire map |
Contributor
There was a problem hiding this comment.
psi_vals is documented above, is this a duplicate?
Comment on lines
+999
to
+1111
| def build_psi_vals_lsq(self, full_currents_vec): | ||
| """ | ||
| Construct a least-squares system enforcing direct flux value constraints. | ||
|
|
||
| This method builds the optimisation system: | ||
|
|
||
| A I_control ≈ b | ||
|
|
||
| where: | ||
|
|
||
| ψ_model(R_i, Z_i) = ψ_target(R_i, Z_i) | ||
|
|
||
| is enforced by matching magnetic flux values at specified locations. | ||
|
|
||
| The optimisation residual is defined as: | ||
|
|
||
| b = ψ_tokamak(I) + ψ_plasma | ||
| - ψ_target | ||
| - ⟨b⟩ | ||
|
|
||
| Mean flux removal is applied to remove arbitrary vertical flux offsets, | ||
| since the Grad–Shafranov equation is invariant under constant flux shifts. | ||
|
|
||
| Parameters | ||
| ---------- | ||
| full_currents_vec : ndarray | ||
| Full coil current vector. | ||
|
|
||
| Example: | ||
| eq.tokamak.getCurrentsVec() | ||
|
|
||
| Returns | ||
| ------- | ||
| A : ndarray | ||
| Jacobian matrix mapping coil current perturbations → flux changes. | ||
|
|
||
| b : ndarray | ||
| Flux mismatch residual vector. | ||
|
|
||
| normalised_loss : list of float | ||
| Normalised constraint violation magnitude. | ||
|
|
||
| Notes | ||
| ----- | ||
| This constraint formulation is commonly used for: | ||
|
|
||
| • Magnetic axis pinning | ||
| • Boundary flux matching | ||
| • Profile shape control | ||
|
|
||
| Mathematical formulation | ||
| ------------------------ | ||
| Solve: | ||
|
|
||
| min_I || G I + ψ_plasma − ψ_target ||² | ||
| """ | ||
|
|
||
| # flux response wrt coil currents | ||
| A_r = self.Gbr_2nd_order[self.control_mask].T | ||
| b_r = np.sum( | ||
| self.Gbr_2nd_order * full_currents_vec[:, np.newaxis], axis=0 | ||
| ) # coils contribution | ||
| b_r += self.Brp_2nd_order # plasma contribution | ||
| loss = [np.linalg.norm(b_r)] | ||
|
|
||
| # Bz field constraint | ||
| A_z = self.Gbz_2nd_order[self.control_mask].T | ||
| b_z = np.sum( | ||
| self.Gbz_2nd_order * full_currents_vec[:, np.newaxis], axis=0 | ||
| ) # coils contribution | ||
| b_z += self.Bzp_2nd_order # plasma contribution | ||
| loss.append(np.linalg.norm(b_z)) | ||
|
|
||
| # dBrdr field constraint | ||
| A_r_deriv = self.Gdbrdr_2nd_order[self.control_mask].T | ||
| b_r_deriv = np.sum( | ||
| self.Gdbrdr_2nd_order * full_currents_vec[:, np.newaxis], axis=0 | ||
| ) # coils contribution | ||
| b_r_deriv += self.dBrdrp_2nd_order # plasma contribution | ||
| loss.append(np.linalg.norm(b_r_deriv)) | ||
|
|
||
| # dBzdz field constraint | ||
| A_z_deriv = self.Gdbzdz_2nd_order[self.control_mask].T | ||
| b_z_deriv = np.sum( | ||
| self.Gdbzdz_2nd_order * full_currents_vec[:, np.newaxis], axis=0 | ||
| ) # coils contribution | ||
| b_z_deriv += self.dBzdzp_2nd_order # plasma contribution | ||
| loss.append(np.linalg.norm(b_z_deriv)) | ||
|
|
||
| # dBrdz field constraint | ||
| A_r_deriv_cross = self.Gdbrdz_2nd_order[self.control_mask].T | ||
| b_r_deriv_cross = np.sum( | ||
| self.Gdbrdz_2nd_order * full_currents_vec[:, np.newaxis], axis=0 | ||
| ) # coils contribution | ||
| b_r_deriv_cross += self.dBrdzp_2nd_order # plasma contribution | ||
| loss.append(np.linalg.norm(b_r_deriv_cross)) | ||
|
|
||
| # dBzdr field constraint | ||
| A_z_deriv_cross = self.Gdbzdr_2nd_order[self.control_mask].T | ||
| b_z_deriv_cross = np.sum( | ||
| self.Gdbzdr_2nd_order * full_currents_vec[:, np.newaxis], axis=0 | ||
| ) # coils contribution | ||
| b_z_deriv_cross += self.dBzdrp_2nd_order # plasma contribution | ||
| loss.append(np.linalg.norm(b_z_deriv_cross)) | ||
|
|
||
| A = np.concatenate( | ||
| (A_r, A_z, A_r_deriv, A_z_deriv, A_r_deriv_cross, A_z_deriv_cross), axis=0 | ||
| ) | ||
| b = -np.concatenate( | ||
| (b_r, b_z, b_r_deriv, b_z_deriv, b_r_deriv_cross, b_z_deriv_cross), axis=0 | ||
| ) | ||
| return A, b, loss | ||
|
|
Contributor
There was a problem hiding this comment.
Has this been added by mistake?
- A function
build_psi_vals_lsqis already defined (just below) - Isn't this a duplicate of
build_null_points_2nd_order_lsq?
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
In addition to being able to set first order magnetic null constraints in the inverse solver, we've now added functionality for setting second order null constraints. This is particularly useful for making snowflake divertor geometries.
There are updates to the Example01b notebook for how to use these constraints and an example snowflake equilibrium.
Note that this requires PR 50 in FreeGS4E.