BUG: Fix flux calculation on surface when there is no fluxing_field

[jzuhone]

PR Summary

Nominally, when one creates a YTSurface object using ds.surface, and then uses calculate_flux on the object to calculate a flux through the surface, if no fluxing_field is provided it should default to 1.0 dimensionless. I discovered by accident today that this isn't the case--mainly because np.ones_like now preserves units (see code below). This is the fix for that. The test has been improved to catch it.

PR Checklist

• New features are documented, with docstrings and narrative docs
• Adds a test for any bugs fixed. Adds tests for new features.

[neutrinoceros]

Looks good to me, just one question: do we know since when np.ones_like started preserving units ? I checked that there is no difference between unyt 2.8 and 2.9.2 (current), using numpy 1.22.3 (current)
There could be a lot of similar bugs, might be interesting to know what changed and when.