Fix logical errors and typos in Spec equilibrium object

[missing-user]

There were a few bugs in the Spec equilibrium object that I spotted:

• Many error checks were written, but the raise clause was omitted, so they didn't take effect
• The RHS of the assignments to si.iota[0:self.nvol+1] and si.oita[0:self.nvol+1] had the wrong shape
• array shape check in as_spec compared tuple to list and always raised
• the initial guess should have nvol elements, not mvol-1

[missing-user]

I noticed, that parts of the codebase assume all profiles to be the same size (mvol), in this example according to the error message local_full_x = self.mvol

simsopt/src/simsopt/mhd/profiles.py

Lines 93 to 98 in 9fc2c36

	# Check that volume index is within bounds
	if (lvol < 0).any():
	raise ValueError('lvol should be larger or equal than zero')
	if (lvol >= self.local_full_x.size).any():
	raise ValueError('lvol should be smaller than Mvol')

or here, where volume_current_profile is simply overwritten in the vacuum region for freeboundary runs, so although nvol elements are used, the profile has to be at least mvol in size:

simsopt/src/simsopt/mhd/spec.py

Lines 1001 to 1019 in 9fc2c36

	if self.volume_current_profile is not None:
	# Volume current is a cumulative profile; special care is required
	# when a dofs is changed in order to keep fixed dofs fixed!
	old_ivolume = copy.copy(si.ivolume)
	for lvol in range(0, self.mvol):
	if self.volume_current_profile.is_fixed(lvol):
	if lvol != 0:
	si.ivolume[lvol] = si.ivolume[lvol] - \
	old_ivolume[lvol - 1] + si.ivolume[lvol - 1]
	self.set_profile(
	'volume_current', lvol=lvol, value=si.ivolume[lvol])
	else:
	si.ivolume[lvol] = self.get_profile('volume_current', lvol)
	if si.lfreebound:
	si.ivolume[self.nvol] = si.ivolume[self.nvol - 1]
	self.volume_current_profile.set(
	key=self.mvol - 1, new_val=si.ivolume[self.nvol - 1])

while other parts create different sizes for different profiles.

simsopt/src/simsopt/mhd/spec.py

Lines 766 to 781 in 9fc2c36

	profile_dict = {
	'pressure': {'specname': 'pressure', 'cumulative': False, 'length': self.nvol},
	'volume_current': {'specname': 'ivolume', 'cumulative': True, 'length': self.nvol},
	'interface_current': {'specname': 'isurf', 'cumulative': False, 'length': self.nvol-1},
	'helicity': {'specname': 'helicity', 'cumulative': False, 'length': self.mvol},
	'iota': {'specname': 'iota', 'cumulative': False, 'length': self.mvol},
	'oita': {'specname': 'oita', 'cumulative': False, 'length': self.mvol},
	'mu': {'specname': 'mu', 'cumulative': False, 'length': self.mvol},
	'pflux': {'specname': 'pflux', 'cumulative': True, 'length': self.mvol},
	'tflux': {'specname': 'tflux', 'cumulative': True, 'length': self.mvol}
	}
	profile_data = self.inputlist.__getattribute__(profile_dict[longname]['specname'])[0:profile_dict[longname]['length']]
	profile = ProfileSpec(profile_data, cumulative=profile_dict[longname]['cumulative'], psi_edge=self.inputlist.phiedge)
	profile.unfix_all()
	self.__setattr__(longname + '_profile', profile)

What is the intended design @smiet ? I think variable size will result in stricter error handling, whereas having all profiles the same size makes them easier to switch out, iterate over, etc. On the other hand, having the same size might confuse users about which elements actually affect the result.

[missing-user]

Also @smiet , shouldn't iota profiles be nvol+1 elements and not mvol? How can the vacuum computational boundary have a prescribed iota?
#457 (comment)

[smiet]

@missing-user Thank you for the great work! Sorry to keep you waiting. The inconsistencies are not intentional and need to be addressed.

The design relies on the user to enforce that the 'ProfileSpec' object is the right size and type for the profile in question. There are some helper functions that aid the user in generating profiles. From the simsopt side it is important that the profiles only contain the degrees-of-freedom that actually affect the plasma, so having them all mvol is not an option.

I noticed, that parts of the codebase assume all profiles to be the same size (mvol), in this example according to the error message local_full_x = self.mvol

The error is in the message here, the test compares to the number of degrees-of-freedom of the SpecProfile instance against the 'volume' (interface label) that is trying to be set, which depends on the profile, and the message should be changed to reflect that.

or here, where volume_current_profile is simply overwritten in the vacuum region for freeboundary runs, so although nvol elements are used, the profile has to be at least mvol in size:

This is due to the (confusing) use of 'cumulative profiles'. The profile on the simsopt side has nvol degrees-of-freedom, but on the SPEC side needs to populate mvol elements of the array.
The value of the volume current in the vacuum needs to be changed if the current inside of it is changed, but this cannot be to an arbitrary value and should not be a DoF.

We have to keep the DoFs on the simsopt side correct, as the optimizer dislikes having degrees-of-freedom that do not change anything (or worse, they do change the elements of the array, leads to calculations that do not reflect physical reality, i.e. current in the vacuum region).

I have not tested the profiles extensively, and there are indeed likely inconsistensies that you point out, like:

Also @smiet , shouldn't iota profiles be nvol+1 elements and not mvol? How can the vacuum computational boundary have a prescribed iota? #457 (comment)

Yes, this must be wrong.

I now see that the the Spec.activate_profile() function might not work correctly for all profiles, volume-current should skip the first element for example.

I see you refactored the NormalField, but there is a simsopt specific issue with this: when changing a single element in the array, the recompute bell is not triggered. Setting the entire Optimizable.local_full_x does trigger it, or you can do it manually.

I did very much appreciate the full-array setters and getters set_vns_vns_asarray as it made code more readable, and it is used in CoilNormalField which gets the full array from a Fourier Transform, and needs it to be set. Could you rename the deleted get_dofs to _array_to_dofs that acts on the input array and sets the local_full_x?

[missing-user]

Thanks for the thorough feedback @smiet , I'll pick this back up after I submit my thesis next month. Just wanted to confirm know that I haven't forgotten about this PR.