Undocumented behavior in f90wrapped SPEC and VMEC (singleton methods)

[smiet]

In discussing with @andrewgiuliani and @mbkumar about the strange behavior of the Spec class (#357), and whether it's fix (#431) could lead to an unsafe memory-leak, it became clear to us that there are certain behaviors of f90wrapped code that are very counter-intuitive and should be documented for the users.

See for example this unexpected behavior:

from simsopt.mhd import Vmec
myvmec1 = Vmec()
# It is initialized with only zeros
print(myvmec1.indata.ac)
>>> [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
# instantiate a new Vmec object
myvmec2 = Vmec()
# modify an element of one of the FORTRAN-side arrays
myvmec2.indata.ac[0]=42
print(myvmec2.indata.ac)
>>> [42.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.]
# BUT! also myvmec1 now has a modified array!
print(myvmec1.indata.ac)
>>> [42.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.]
# also the base class, that we now import already has the modified array
import vmec
print(vmec.vmec_input.ac)
>>> [42.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.]

Changing the ac attribute in myvmec2 has also changed that attribute in myvmec1 (and in the base FortranModule vmec.vmec_input).

This is true for any f90wrapped FortranModule instance, as they use a __metaclass__ = Singleton- type construction that returns the existing class if called more than once.

This behavior seems pretty foundational to the design of f90wrap (the code implementing it was committed 10 years ago, at the very start), though it is not well-documented why this was chosen, and took me by surprise initially. This is also one of the reasons that the CI was failing; each test instantiating a new Spec() object basically gets returned the already existing class.

If only one class is active at the same time, I don't see an immediate issue as both SPEC and VMEC (and FORTRAN in general) allocate and define all arrays during program execution, so on the FORTRAN side every used array will be correct. If more than one object is created on the python side, things can go very wrong, as they start writing each others' memory, and I have seen segmentation faults and crashes occur.

We should also investigate if this behavior can lead to memory leaks, and if there could be any issues from this behavior when running in parallel, both through openMP and MPI

But primarily this should be documented on our side for both the Vmec and Spec class. In fact, I see the documentation contains no top-level entries on the MHD codes:

I am unfamiliar in how to create these structures in sphinx, @mbkumar or @landreman, could you create an MHD entry with sub-entries VMEC200 and SPEC (and possibly DESC, but that is wip?) I will then populate the SPEC description.

I am thinking of a warning of the type:

.. admonition:: warning
   :title: Warning

   Due to a limitation of the fortran-bindings, there can be only one :obj:`~simsopt.mhd.Spec` object active at the same time. Code like
    >>> from simsopt.mhd import Spec
    >>> myspec1 = Spec()
    >>> myspec2 = Spec('my_input.sp')
  will lead to errors and crashes.