Skip to content
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.

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

Planar coil array branch merge #464

Open
wants to merge 202 commits into
base: master
Choose a base branch
from
Open

Planar coil array branch merge #464

wants to merge 202 commits into from

Conversation

akaptano
Copy link
Contributor

Hi all,

Here is a branch I have that is built off of some of the machinery from Siena's coil forces branch. Some notes:

  1. The branch implements various pointwise and net force and torque calculations, and the total vacuum energy.
  2. The branch implements a JaxCurvePlanarFourier class and fixes the Jacobian calculation (thanks Alex) for the CurvePlanarFourier class, which had a bug.
  3. It adds some geometry-related functions and various other functionality for effectively optimizing arrays of dipole coils, including with forces and torques.
  4. There is still a fair bit of stuff in the 3_Advanced/ examples directory to reproduce the paper I'm putting on arXiv shortly, which I can put its own Zenodo repo if need be. I also included Siena's Pareto-scan scripts, since I actually think future coil optimizers would find these helpful to not have to recreate.
  5. There are some lingering issues. As the number of coils gets large (~40), the number of Optimizable weak references seem to spiral out of control in any of the force/torque objectives that need to generate their own BiotSavart objects. I have been able to somewhat control this by regularly erasing weak references when making calls to this terms. This seems to have no bearing on the optimization itself. I have also not yet correctly implemented the fast version of the self inductance for a rectangular coil from Siena and Matt.
  6. I am planning on adding passive coil functionality, but figure I should stop here and continue on a branch split off from this one.

Let me know if you have any questions, and thanks for any time you can spare to take a look!

…arted) geometry file for the grid of coils. Next step is to fill this out and write out a list of coils to curve files which can be visualized in Paraview.
…nloading simsopt with recursive and hopefully the c++ backend installation will work better.
…on troubles) into a geometry class in python. self and mutual inductance calculations seem to be working well on simple unit tests. Wrote a script for unit testing the inductances.
…tion disagree for two-coil checks with a rotation alpha = 0.5pi and delta = 0, but inductance looks correct, since this should be equivalent to the analytic formula with coils offset but with the same axis of orientation, so somehow the flux calculation must be wrong. Not sure how... I have now exhaustively checked that the integration points in the integral are correct, that B is correct, that the normal vector n to the plane is correct, etc.
… are working just fine -- the issue was the inductances calculation. I need to account for both (1) the separate rotations of coils i and coils j, and (2) the center of coil i should be subtracted off from X, Y, Z in the calculation, since the origin is supposed to be at this center.
…s are visualized on the coils. Inductance calculation was slow so wrote a script to compute it in c++, and got the c++ compiling correctly. Need to make sure the calculation is correct since tests are not yet passing. Added calculation for the magnetic fields and currents from all the coils.
…ch faster than the python version. However, not clear if this can be used in optimization. Perhaps for auto differentiation of L_inv we can use the python version, but for direct evaluations of L we can use the c++ version. Unclear if Jax will be okay with the python-wrapped version of the c++ function.
…ecause all the work is setting the points in the call to BiotSavart, so better to do this in python. Moreover, the convergence of the integral is slow, so using dblquad from scipy produces much better results than doing naive uniform or gauss-legendre integration.
…l to calculate the fluxes is wayyyy to slow, so need to speed this up. Tried writing up a version of this in c++ but realized I actually calculated the fluxes from the PSCs, which is not what is needed. We need to have a fast calculation of the fluxes from the fixed TF coils.
…s for the flux calculation. Plan is to then set the Bfield as usual, then call another c++ function to do the integration.
…rix in optimization. This is tricky since I think there is a double Rodrigues rotation around the normal vector of each coil in the calculation. Agrees now with the calculation from the existing CircularCoil class, but for some reason the pure BiotSavart calculation disagrees. Need to square away all these rotations and put this issue to bed + add a unit test to check the three Bn calculations. Then, will run a basic optimization again and make a little visual. Then next step will be to derive and code up the least squares jacobian.
…tion appears to be very simple. The Euler rotation matrix simply is applied to B so that under rotation we get RB(RT * x). It doesnt matter if this matrix is the Euler one or the fancier one used in the CircularCoil code. Should add comments to clarify this for future users. Code is an absolute mess right now because of all my debugging but thought it was worth saving in this state for posterity. Next step is to try and deal with an issue that the inductance matrix appears very ill-conditioned so gives dumb solutions for the currents. Not sure if this is because there is a bug in the calculation or inherent in the problem.
…e I had initialized a stupid grid where a few of the coils were touching or intersecting, so that indeed the mutual inductances were off the charts. After initializing a grid where the coils are well-spaced, this problem evaporates and in fact L is very well-conditioned since it is mostly dominated by the self-inductances.
… Bn in c++ compared with direct biotsavart and circularCoil classes. Inductance and Bn calculations are < 1 second now with almost 3000 PSCs but still slower than I would like. Nonetheless, next step appears to be coding up the Jacobians.
…code. Seems pretty good so far, but cant seem to figure out encoding the reflection symmetry in the quaternions, which are typically used for proper rotations. This is anyways probably not worth it, since can always call coils_via_symmetries function or apply_symmetries_to_curves function. Saving now before doing some cleanup.
…now that stellarator symmetry and field period symmetry is imposed correctly when calculating the magnetic fields. So far not quite getting the right answer, but all the machinery is there. Need to convince myself about what symmetries do to the inductance matrix, since I am not sure about how the current matrix transforms. Usually the currents in the coils are imposed, so you can simply take a set of N coils and apply the symmetries to get the rest. Here, the currents in the coils are themselves dependent on the coils that you get after you symmetrize.
…at is going on with the stellarator and field period symmetries. Made some intuition progress but lots of work still to do. Not sure if the fluxes should be equal in the symmetrized coils but direct calculation appears to suggest that.
…ely. However, looks like the Bn calculation from the c++ code does not agree with the same calculation from direct BiotSavart when there are discrete symmetries, so still need to debug this before I can conclude that I am actually doing the optimization correctly for QH.
…ries. Still a few sign discrepancies and lots of clean up to do.
… fiddling to get them looking more plausible but still no luck. Going to do some code cleanup to take a little break.
… in one coil configuration. Just testing the A matrix gradient so far and havent been able to quite get it working, even though I can reproduce the A matrix itself. Need to think harder about what is returned when I call dB_by_dX.
…. Serious addressing of the derivative calculations will need to wait for now.
…g. Had to fix alphas and deltas to be between -pi and pi to play nice with arctan2 function in numpy and consistently avoid changing the angles when I manipulate them. Started writing up explicit function to compute dB_dkappa which is needed to compute derivatives of A and psi part of objective.
…the passive coils. Already have dL_dalpha so hopefully will finish this up tomorrow and then start on the derivative for psi. Added unit test for dA_dalpha, and wrote another function A_matrix_explicit to double check some of the rotations going on in the explicit biot savart calculation. A_matrix_explicit agrees well with A_matrix calculation, but dA_dalpha still not quite right.
…alled with the wrong argument in the python code. Now trying to run all the unit tests and seeing some sign errors with respect to Bn computed with the CircularCoil class. My guess is that sign errors from the arctan2 calculations are causing this. Need to fix the PSC grid class to not keep converting between normal vectors and the alpha and delta maybe.
…o getting all the unit tests correct is the issue converting between the normal vector and the Euler angles.
@akaptano akaptano self-assigned this Dec 16, 2024
@akaptano akaptano added the enhancement New feature or request label Dec 16, 2024
akaptano and others added 6 commits December 17, 2024 09:45
…the pareto runs. Checked that it looks like tve is being optimized reasonably. Going to try another set of small runs to verify if anything changes.
…being used without the signs in LijIiIj. Rerunning the TVE scan.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

1 participant