R_modes = np.vstack(
    (
        [0, 0, 0],
        eq.surface.R_basis.modes[
            np.max(np.abs(eq.surface.R_basis.modes), 1) > k, :
        ],
    )
)
Z_modes = eq.surface.Z_basis.modes[
    np.max(np.abs(eq.surface.Z_basis.modes), 1) > k, :
]
boundary constraints = (
    FixBoundaryR(eq=eq, modes=R_modes),
    FixBoundaryZ(eq=eq, modes=Z_modes),
)

def multigrid_optimize(things, objective, constraints, ms):
    for m in ms:
        cons = constraints = (FixBoundaryR(modes > m), FixBoundaryR(modes > m))
        optimizer.optimize(things, objective, cons)