adapt BS pointsetting in QFM label functions

[smiet]

QfmSurface can take a label, which is an Optimizable that constrains the surface to a desired quantity. Often used are Volume and ToroidalFlux.

There is an issue that surfaceobjectives ToroidalFlux and QfmResidual set the BiotSavarts evaluation points to differently shaped arrays; QfmResidual needs to evaluate on the entire surface, whereas ToroidalFlux only on a poloidal loop.
The re-setting of the evaluation points was only done when the recompute_bell was rung, allowing re-use of intermediate results calculated in the BiotSavart object, if for example derivatives are requested.

But when optimizing QFM surfaces for Toroidal flux, the points would still be set to the full surface from the b.n calculation, causing errors.

In the test_qfm.py we used to create a separate bs_tf object for the ToroidalFlux evaluation to avoid this, but this was confusing: a single BiotSavart should be able to provide all the field-evaluation needs.

fix:

The surfaceobjectives now have an attribute need_to_invalidate, which is set when the recompute bell is rung.
They also have an attribute correct_shape which is the shape of points the BiotSavart expects.
At each evaluation the shape is tested against correct_shape, re-setting the points only if the shape is not correct.

why?

This should maintain the speed, as the cache is not invalidated at every evaluation (Evaluations of the derivatives of A and B can reuse previous computations).

MWE of previously-failing code:

from simsopt.geo import SurfaceRZFourier
from simsopt.geo.qfmsurface import QfmSurface
from simsopt.geo.surfaceobjectives import QfmResidual, ToroidalFlux
from simsopt.configs import get_data


base_curves, base_currents, ma, nfp, bs = get_data('w7x')

s = SurfaceRZFourier.from_nphi_ntheta(mpol=8, ntor=8, stellsym=True, nfp=nfp, range='half period')

s.fit_to_curve(ma, 0.15, flip_theta=True)

label=ToroidalFlux(s, bs)
toroidal_flux_target=20
qfm_surface = QfmSurface(bs, s, label, toroidal_flux_target)

qfm_surface.minimize_qfm_penalty_constraints_LBFGS(tol=1e-20, maxiter=50,
                                                   constraint_weight=1e-3)

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 89.99%. Comparing base (b1eb111) to head (81bf2cc).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #594      +/-   ##
==========================================
+ Coverage   89.97%   89.99%   +0.02%     
==========================================
  Files          84       84              
  Lines       17364    17382      +18     
==========================================
+ Hits        15623    15643      +20     
+ Misses       1741     1739       -2     

Flag	Coverage Δ
unittests	89.99% <100.00%> (+0.02%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[smiet]

should be a very simple PR, @akaptano or @mishapadidar, can you do the honors of approving?

[akaptano]

@smiet can you address the comment I left? I will be happy to approve after.

[andrewgiuliani]

hi @smiet, I am confused, why wasn't

    def recompute_bell(self, parent=None):
        self.invalidate_cache()

sufficient?

apparently the optimizable depends on biotsavart so if the points in the biotsavart object have been changed, the then the recompute_bell should invalidate the cache

[smiet]

'BiotSavart' does not depend on its points, only on the dofs of coils. Should it?

[andrewgiuliani]

simsopt/src/simsopt/field/magneticfield.py

Line 27 in b1eb111

The cache is automatically cleared when ``set_points`` is called or one of the dependencies

the cache should be cleared when set_points() is called, which I would assume calls recompute_bell then invalidate_cache() of the surfaceobjective

[smiet]

That's maybe The issue, set_points() does not call the recompute bell. I also do not think it should. Evaluating the field at a new point does not warrant waking up the whole continent. The BS object is happy to keep providing it's services to whomever wants to evaluate a field wherever, nothing fundamentally had changed.

[andrewgiuliani]

but isn't this the issue? you want the surfaceobjective cache to be cleared when the points in the biotsavart object are set, so that when the surfaceobjective needs the B field, it will know to recompute everything at the right points.

[andrewgiuliani]

that's why I'm confused:

simsopt/src/simsopt/field/magneticfield.py

Line 27 in b1eb111

The cache is automatically cleared when ``set_points`` is called or one of the dependencies

when set_points() is called, everything should be cleared, so I think that means we have to understand why the current implementation doesnt work. I don't see why it's not working, and it'll take some digging.

[andrewgiuliani]

My confusion is: set_points should automatically be calling the recompute bell, but clearly it's not. Hope I'm explaining

[smiet]

The C++ side cache is cleared, yes, but it does not state that the dependency graph on the python side is called into action. Many functions on the python side rely on changing the evaluation points over and over, which does not alter any fundamental property of the underlying mathematical Optimizable

[smiet]

I do not know if it's desired that set_points always calls the recompute_bell(). Set_points is called thousands of times in for example integrators, and this would add a lot of overhead. I think we should expose our blindingly fast Biot-Savart without any hindrances to any calling function, but I'm biased, because pyoculus relies on this.

[andrewgiuliani]

To discuss at the next dev meeting so we can chat 👍🏻

[smiet]

Sorry, @akaptano I cannot see any comment of yours, it's it still 'pending'?

[akaptano]

Sorry to ask, but could you add good docstrings to this file (using Cursor is fine) while you're at it? For instance, I am not sure about the line self.correct_shape = self.surface.gamma()[self.idx].shape since I don't know what idx does. And I thought that self.surface.gamma() should basically always shape (nphi, ntheta, 3), so isn't self.surface.gamma()[0, :, :].shape the same as self.surface.gamma()[1, :, :].shape, etc.?