Cannot load and re-solve old DESC files with the current master.

[rahulgaur104]

I was trying to load an old NCSX DESC file and re-solve it with DESC using the script below:

from desc import set_device
set_device("gpu")

import numpy as np
from desc.equilibrium import EquilibriaFamily, Equilibrium
from desc.grid import LinearGrid, Grid
from desc.objectives import (
    ForceBalance,
    ObjectiveFunction,
)
from desc.vmec import VMECIO
from desc.backend import jnp
from desc.equilibrium.coords import compute_theta_coords

fam = EquilibriaFamily()

eq = EquilibriaFamily.load("NCSX_og1_reinit.h5")[-1]

eq.L = 16
eq.M = 16
eq.N = 16

eq, result = eq.solve(
    objective="force",
    ftol=5e-4,
    xtol=1e-6,
    gtol=1e-6,
    maxiter=10,
    verbose=3,
    copy=True,
)

which throws the following error

/scratch/gpfs/rg6256/DESC_fork/DESC/desc/io/hdf5_io.py:125: RuntimeWarning: Save attribute '_rho' was not loaded.
  warnings.warn(
/scratch/gpfs/rg6256/DESC_fork/DESC/desc/io/hdf5_io.py:125: RuntimeWarning: Save attribute '_shift' was not loaded.
  warnings.warn(
/scratch/gpfs/rg6256/DESC_fork/DESC/desc/io/hdf5_io.py:125: RuntimeWarning: Save attribute '_rotmat' was not loaded.
  warnings.warn(

Building objective: force
Precomputing transforms
Timer: Precomputing transforms = 3.14 sec
Timer: Objective build = 13.3 sec

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[2], line 28
     25 eq.M = 16
     26 eq.N = 16
---> 28 eq, result = eq.solve(
     29     objective="force",
     30     ftol=5e-4,
     31     xtol=1e-6,
     32     gtol=1e-6,
     33     maxiter=10,
     34     verbose=3,
     35     copy=True,
     36 )
     38 #eq = EquilibriaFamily.load("eq_final.h5")
     39 #print(len(eq))

File /scratch/gpfs/rg6256/DESC_fork/DESC/desc/equilibrium/equilibrium.py:1794, in Equilibrium.solve(self, objective, constraints, optimizer, ftol, xtol, gtol, maxiter, x_scale, options, verbose, copy)
   1788 if self.bdry_mode == "poincare":
   1789     raise NotImplementedError(
   1790         "Solving equilibrium with poincare XS as BC is not supported yet "
   1791         + "on master branch."
   1792     )
-> 1794 things, result = optimizer.optimize(
   1795     self,
   1796     objective,
   1797     constraints,
   1798     ftol=ftol,
   1799     xtol=xtol,
   1800     gtol=gtol,
   1801     x_scale=x_scale,
   1802     verbose=verbose,
   1803     maxiter=maxiter,
   1804     options=options,
   1805     copy=copy,
   1806 )
   1808 return things[0], result

File /scratch/gpfs/rg6256/DESC_fork/DESC/desc/optimize/optimizer.py:192, in Optimizer.optimize(self, things, objective, constraints, ftol, xtol, gtol, ctol, x_scale, verbose, maxiter, options, copy)
    190 if len(linear_constraints):
    191     objective = LinearConstraintProjection(objective, linear_constraints)
--> 192     objective.build(verbose=verbose)
    193     if nonlinear_constraint is not None:
    194         nonlinear_constraint = LinearConstraintProjection(
    195             nonlinear_constraint, linear_constraints
    196         )

File /scratch/gpfs/rg6256/DESC_fork/DESC/desc/optimize/_constraint_wrappers.py:82, in LinearConstraintProjection.build(self, use_jit, verbose)
     80 for con in self._constraints:
     81     if not con.built:
---> 82         con.build(verbose=verbose)
     84 self._dim_f = self._objective.dim_f
     85 self._scalar = self._objective.scalar

File /scratch/gpfs/rg6256/DESC_fork/DESC/desc/objectives/linear_objectives.py:273, in BoundaryRSelfConsistency.build(self, use_jit, verbose)
    270 if eq.bdry_mode == "lcfs":
    271     j = np.argwhere((modes[:, 1:] == [m, n]).all(axis=1))
    272     surf = (
--> 273         eq.surface.rho
    274         if self._surface_label is None
    275         else self._surface_label
    276     )
    277     self._A[j, i] = zernike_radial(surf, l, m)
    278 else:

File /scratch/gpfs/rg6256/DESC_fork/DESC/desc/geometry/surface.py:161, in FourierRZToroidalSurface.rho(self)
    158 @property
    159 def rho(self):
    160     """float: Flux surface label."""
--> 161     return self._rho

AttributeError: 'FourierRZToroidalSurface' object has no attribute '_rho'

Is anyone else seeing the same error? Aren't there any tests to check this?

[ddudt]

I think you just have to set eq.surface._rho = 1 after you load the equilibrium, and then it should run fine. This is because of the changes made in PR #814

We do not officially guaranteed the code to be backwards compatible, but usually it's only minor changes like this that are easy to fix manually.

[rahulgaur104]

I understand the definition of backwards compatibility. However, for a non-developer, it is by no means clear how to get rid of that error. Please make sure to add a warning, or a test that can be turned on whenever v1.0 is released.