New refit without contouring

[bclyons12]

This change bypasses the standard MXH fitting infrastructure

Old steps

1. Map Psi to a uniform (R, Z) grid
2. Contour the surfaces to get (R, Z) points at each desired Psi
3. Use an algorithm very similar if not identical to GACODE to find the coefficients that match the set of contour points

New steps

1. Find the minimum and maximum (R, Z) points that match the desired Psi (using JuMP to constrain the Psi value well)
2. Knowing Zmin and Zmax, I know the Z values that correspond to equally spaced theta
3. Use Roots to find the R values at each desired Z that has the desired Psi value
4. Compute the Theta_r in the MXH R definition for each point
5. Since this a uniform grid in theta, use any FFT to find the Fourier coefficients in Theta_r

[orso82]

Does this make use of the recent auto differentiation capabilities in TEQUILA?

[jcandy]

Is there a comparison yet of profile variation of shape coefficients between GACODE and the various FUSE components? I would like to see this.

[bclyons12]

@orso82 Yes! I just merged #35 into master allowing this

[bclyons12]

@jcandy Not yet, but we could do it. Do you have a specific case you'd like to look at?

[orso82]

I have used JuMP in the past and it worked.

However I find the metaprogramming associated with it to be a bit clunky. I am thinking that perhaps we could easily write a optimization function that has a similar API to Optim.jl but that runs JuMP underneath. Could this be a good way to systematically move away from using Optim?

This thought aside, looks good to me and I have no comments. Merge when you are ready.

[bclyons12]

@jcandy See the edited description at the top of the issue now that describes what's happening here.

[orso82]

@bclyons12 maybe comment why you chose to go with this new method ;)

[bclyons12]

@orso82 Yes, JuMP is a little weird, but I found it's great for nonlinear constraints. I need pretty much the exact (R,Z) extrema points and JuMP gives it to me. I'm not sure JuMP is optimal, but this is so fast [milliseconds instead of O(second) before] I don't think we need to worry about it.

I'm just testing for robustness now, but so far, so good.

[bclyons12]

@orso82 Good call.

The old method is fundamentally limited by its Step 1. You have to put Psi onto a grid, so the accuracy of your flux surfaces is limited by the resolution of the grid. Calls to Psi(R,Z) are somewhat slow in TEQUILA, since you need to do the nonlinear conversion of (R, Z) to (rho, theta), so going to high resolution in the grid was dominating the runtime in FUSE. Going to lower resolution sped things up, but limits the amount of convergence you can get.

This takes care of all of that. find_extrema operates in (rho, theta) space, so it doesn't need to do the nonlinear conversion from (R, Z). Then I find the minimum number of points that "exactly" match the desired Psi on the contour. This should let me arbitrarily converge the accuracy of flux surfaces while being quite fast.

Here are results running at moderate resolution for

shot =  Shot(21, 21, 15, "sample/g_chease_mxh_d3d");
@time refill = solve(shot, 50; debug=true);

Old

Iteration 1: Ψaxis = -2.2858782969834026, Error: 0.41487725858265145
Iteration 2: Ψaxis = -2.0317544999432084, Error: 0.25412379704019417
Iteration 3: Ψaxis = -1.9351054911146859, Error: 0.09664900882852256
Iteration 4: Ψaxis = -1.8963719064589457, Error: 0.03873358465574017
Iteration 5: Ψaxis = -1.8811899359602058, Error: 0.015181970498739927
Iteration 6: Ψaxis = -1.875202356864803, Error: 0.005987579095402706
Iteration 7: Ψaxis = -1.8729428023539025, Error: 0.002259554510900541
Iteration 8: Ψaxis = -1.871606931709013, Error: 0.001335870644889603
Iteration 9: Ψaxis = -1.8710183906576507, Error: 0.0005885410513621903
Iteration 10: Ψaxis = -1.8709404766420774, Error: 7.791401557333266e-5
Iteration 11: Ψaxis = -1.870872642285798, Error: 6.783435627943923e-5
Iteration 12: Ψaxis = -1.870828160076175, Error: 4.448220962305349e-5
Iteration 13: Ψaxis = -1.8708169552885054, Error: 1.1204787669472083e-5
Iteration 14: Ψaxis = -1.8708108613065368, Error: 6.0939819686023355e-6
Iteration 15: Ψaxis = -1.8708120320867423, Error: 1.1707802054505834e-6
Iteration 16: Ψaxis = -1.8708196832958777, Error: 7.651209135373094e-6
Iteration 17: Ψaxis = -1.8708475939249494, Error: 2.7910629071703497e-5
Iteration 18: Ψaxis = -1.870926091313032, Error: 7.849738808252127e-5
Iteration 19: Ψaxis = -1.870710699131431, Error: 0.0002153921816008264
Iteration 20: Ψaxis = -1.8705909788382367, Error: 0.00011972029319440303
Iteration 21: Ψaxis = -1.870819890687945, Error: 0.00022891184970830203
Iteration 22: Ψaxis = -1.870792341849087, Error: 2.7548838857915925e-5
Iteration 23: Ψaxis = -1.8706993602632387, Error: 9.29815858483618e-5
Iteration 24: Ψaxis = -1.8706791585519087, Error: 2.0201711329992378e-5
Iteration 25: Ψaxis = -1.870737267858489, Error: 5.810930658034508e-5
Iteration 26: Ψaxis = -1.870754861250102, Error: 1.7593391612891196e-5
Iteration 27: Ψaxis = -1.8708174452612376, Error: 6.258401113568013e-5
Iteration 28: Ψaxis = -1.8717928307688738, Error: 0.0009753855076362061
Iteration 29: Ψaxis = -1.8707807641230718, Error: 0.0010120666458020011
Iteration 30: Ψaxis = -1.870092816800134, Error: 0.0006879473229377542
Iteration 31: Ψaxis = -1.8706978149842906, Error: 0.0006049981841564911
Iteration 32: Ψaxis = -1.87065699156786, Error: 4.082341643063536e-5
Iteration 33: Ψaxis = -1.8705873362645424, Error: 6.965530331748404e-5
Iteration 34: Ψaxis = -1.8709203818989315, Error: 0.0003330456343890731
Iteration 35: Ψaxis = -1.871259780887274, Error: 0.0003393989883424542
Iteration 36: Ψaxis = -1.870589967704447, Error: 0.0006698131828268838
Iteration 37: Ψaxis = -1.8703740023039057, Error: 0.00021596540054136248
Iteration 38: Ψaxis = -1.8708959933881626, Error: 0.0005219910842568787
Iteration 39: Ψaxis = -1.8707426337769344, Error: 0.00015335961122819874
Iteration 40: Ψaxis = -1.870582929125762, Error: 0.00015970465117232102
Iteration 41: Ψaxis = -1.8711386218362496, Error: 0.0005556927104874987
Iteration 42: Ψaxis = -1.8707163907167503, Error: 0.00042223111949923897
Iteration 43: Ψaxis = -1.8704771153920552, Error: 0.0002392753246951429
Iteration 44: Ψaxis = -1.870800830686176, Error: 0.0003237152941208965
Iteration 45: Ψaxis = -1.870746443445805, Error: 5.438724037110099e-5
Iteration 46: Ψaxis = -1.8707202139490475, Error: 2.6229496757457227e-5
Iteration 47: Ψaxis = -1.8709868311084903, Error: 0.0002666171594427613
Iteration 48: Ψaxis = -1.8707603352434212, Error: 0.00022649586506906516
Iteration 49: Ψaxis = -1.8705770073087122, Error: 0.00018332793470898778
Iteration 50: Ψaxis = -1.8709735974870432, Error: 0.00039659017833093557
 59.228807 seconds (1.55 M allocations: 2.244 GiB, 0.23% gc time)

New

Iteration 1: Ψaxis = -2.2858782969834026, Error: 0.41487725858265145
Iteration 2: Ψaxis = -2.03111762738907, Error: 0.2547606695943325
Iteration 3: Ψaxis = -1.934677180203255, Error: 0.09644044718581513
Iteration 4: Ψaxis = -1.8961856278945972, Error: 0.03849155230865775
Iteration 5: Ψaxis = -1.8810263750558216, Error: 0.015159252838775616
Iteration 6: Ψaxis = -1.874999577510467, Error: 0.006026797545354512
Iteration 7: Ψaxis = -1.8725551220751664, Error: 0.0024444554353006964
Iteration 8: Ψaxis = -1.871530000287092, Error: 0.0010251217880743457
Iteration 9: Ψaxis = -1.8710822101113864, Error: 0.0004477901757056202
Iteration 10: Ψaxis = -1.8708774345560049, Error: 0.00020477555538156444
Iteration 11: Ψaxis = -1.8707791614580327, Error: 9.827309797216799e-5
Iteration 12: Ψaxis = -1.8707296604043413, Error: 4.9501053691392016e-5
Iteration 13: Ψaxis = -1.8707035931069549, Error: 2.606729738641178e-5
Iteration 14: Ψaxis = -1.8706893099693502, Error: 1.4283137604653717e-5
Iteration 15: Ψaxis = -1.8706812010194303, Error: 8.108949919982678e-6
Iteration 16: Ψaxis = -1.8706764497217079, Error: 4.751297722371817e-6
Iteration 17: Ψaxis = -1.8706735948125186, Error: 2.8549091892760003e-6
Iteration 18: Ψaxis = -1.8706718446600574, Error: 1.7501524611773789e-6
Iteration 19: Ψaxis = -1.8706707514381258, Error: 1.0932219316472924e-6
Iteration 20: Ψaxis = -1.870670065273656, Error: 6.861644696876112e-7
Iteration 21: Ψaxis = -1.8706696284638507, Error: 4.3680980543747694e-7
Iteration 22: Ψaxis = -1.8706693466768733, Error: 2.817869773075188e-7
Iteration 23: Ψaxis = -1.8706691651475733, Error: 1.815293000362317e-7
Iteration 24: Ψaxis = -1.8706690476952814, Error: 1.1745229189230599e-7
Iteration 25: Ψaxis = -1.8706689714747573, Error: 7.622052411448976e-8
Iteration 26: Ψaxis = -1.8706689218807424, Error: 4.959401489479376e-8
Iteration 27: Ψaxis = -1.8706688895635815, Error: 3.231716094731496e-8
Iteration 28: Ψaxis = -1.870668868840324, Error: 2.0723257554422503e-8
Iteration 29: Ψaxis = -1.870668854255342, Error: 1.4584981888887683e-8
Iteration 30: Ψaxis = -1.870668845971885, Error: 8.283457031410535e-9
Iteration 31: Ψaxis = -1.8706688399742726, Error: 5.997612406716257e-9
Iteration 32: Ψaxis = -1.8706688360292574, Error: 3.9450152033992936e-9
Iteration 33: Ψaxis = -1.8706688335495054, Error: 2.479751959327814e-9
Iteration 34: Ψaxis = -1.8706688318978784, Error: 1.6516270573418979e-9
Iteration 35: Ψaxis = -1.8706688308368147, Error: 1.0610636813623842e-9
Iteration 36: Ψaxis = -1.8706688301385337, Error: 6.982809885869301e-10
Iteration 37: Ψaxis = -1.8706688296926988, Error: 4.458349245339832e-10
Iteration 38: Ψaxis = -1.8706688293526084, Error: 3.4009040028593063e-10
Iteration 39: Ψaxis = -1.8706688291763107, Error: 1.762976431507468e-10
Iteration 40: Ψaxis = -1.8706688290407352, Error: 1.3557555078591577e-10
Iteration 41: Ψaxis = -1.8706688289648707, Error: 7.586442585250097e-11
Iteration 42: Ψaxis = -1.8706688289122224, Error: 5.26483301399594e-11
Iteration 43: Ψaxis = -1.8706688288767186, Error: 3.550382210448788e-11
Iteration 44: Ψaxis = -1.8706688288533972, Error: 2.3321344855276038e-11
Iteration 45: Ψaxis = -1.8706688288380002, Error: 1.539701699471152e-11
Iteration 46: Ψaxis = -1.8706688288279771, Error: 1.0023093466315913e-11
Iteration 47: Ψaxis = -1.8706688288213587, Error: 6.618483539000408e-12
Iteration 48: Ψaxis = -1.8706688288170068, Error: 4.351852211925689e-12
Iteration 49: Ψaxis = -1.8706688288142261, Error: 2.780664587476167e-12
Iteration 50: Ψaxis = -1.870668828812366, Error: 1.8600676554569873e-12
 31.436260 seconds (4.20 M allocations: 2.255 GiB, 0.58% gc time)

[orso82]

Fantastic! Could the accuracy be an input parameter. It looks like now you can converge to 1E-6 with much fewer iterations. That should speedup things!

[jcandy]

@jcandy Not yet, but we could do it. Do you have a specific case you'd like to look at?

I sent email on Feb 27 with results for an NSTX case that you sent me

[lstagner]

Will this new approach be used in ExtendedMillerHarmonic.jl?

[bclyons12]

@lstagner This is a different capability from what's in MillerExtendedHarmonic.jl, which has functions for reading in a bunch of (R, Z) points and giving an MXH fit to those points. This method is for getting MXH coefficients such that the curve satisfies Psi(R,Z) = Psi0. We could generalize that to any F(R,Z) = F0 and move it into MillerExtendedHarmonic.jl, creating an MXH contouring capability.

[orso82]

This sounds like a good idea @bclyons12 !

[orso82]

And great point @lstagner !!

[lstagner]

I also think the fft method to determine the coefficients (instead of the current integral method) may be more faster/stable

[bclyons12]

@lstagner agreed when contouring from a field, but it won't work if you just give (R,Z) coordinates, since you need equal spacing in theta to do an FFT

[lstagner]

My thought was to create a spline and then do equal spacing. I do wonder what the spacing would need to be to resolve higher modes. Does the moment method work better for calculating the coefficients?

[bclyons12]

My guess is that the moment method will give the best fit to arbitrarily spaced points. I'm not sure there would be an advantage to splining, downsampling, then taking an FFT. Perhaps it would be faster but that's not obvious to me. Fitting a contour to a field, like here, is a different problem.