Support loading CGYRO simulations run off input.gacode

[cholland]

Hi all- apologies if this feature is supported but I missed it. I think it would be very useful if pyrokinetics supported importing CGYRO simulations run with using PROFILE_MODEL=2 (i.e. read in profiles from an input.gacode file and determine gradients from there, rather than explicitly specifying in the input.cgyro file with PROFILE_MODEL=1). Insofar as the pygacode tools support reading in and plotting simulations run this way I think it should be a pretty doable addition. This capability is something I'd like to utilize in the SMARTS SciDAC project, happy to work with team to implement.

[bpatel2107]

We sort of support this as pyrokinetics can read in input.gacode files. This means we can populate our Equilibrium and Kinetics objects. However pyrokinetics will determine the gradients and local flux surface values directly. This should in theory be similar to the pygacode value, but likely will lead to subtly different results due to interpolation issues.

For example the following would read the input.gacode file and load the local values at psi_n=0.5

# Equilibrium file
eq_file = "input.gacode"

# Kinetics data file
kinetics_file = "input.gacode"

# Load up pyro object
pyro = Pyro(
    eq_file=eq_file,
    eq_type="GACODE",
    kinetics_file=kinetics_file,
    kinetics_type="GACODE",
)

# Generate local parameters at psi_n=0.5
pyro.load_local(psi_n=0.5, local_geometry="MXH")

The issue I think would be when when PROFILE_MODEL= 2 is set as we don't technically read in the input.gacode file though it should be straightforward to read this in and generate the full profile and such. The issue is whether or not to use the pygacode values for the inputs or the pyrokinetics values (we ideally want to match the simulation so we ought to use what CGYRO actually used.). If we want to use the pygacode values then we will need to find a way to get those values out of the CGYRO and into the pyro object. This should be possible using the out.cgyro.info file

So I guess it depends on what you want really. If you are fine to use the code as above then we don't need to do anything. If you want to set PROFILE_METHOD=2 and use the internal interpolation then we need to tweak how we get the physics parameters. I don't see this being a massive issue though.

[jmcclena]

Thanks @bpatel2107. @cholland We need to test how similar the calculated gradients are between pyrokinetics and cgyro. I personally worry that they might differ in the pedestal region, depending on how the gradients are calculated.

@bpatel2107 For you script, I'm getting the follow error on a test input.gacode file:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[3], line 13
     10 kinetics_file = base_dir+"input.gacode"
     12 # Load up pyro object
---> 13 pyro = Pyro(
     14     eq_file=eq_file,
     15     eq_type="GACODE",
     16     kinetics_file=kinetics_file,
     17     kinetics_type="GACODE",
     18 )
     20 # Generate local parameters at psi_n=0.5
     21 pyro.load_local(psi_n=0.5, local_geometry="MXH")

File [~/programming/python/pyrokinetics/src/pyrokinetics/pyro.py:198](http://localhost:8889/lab/tree/~/programming/python/pyrokinetics/src/pyrokinetics/pyro.py#line=197), in Pyro.__init__(self, eq_file, eq_type, eq_kwargs, kinetics_file, kinetics_type, kinetics_kwargs, gk_file, gk_output_file, gk_code, gk_type, nocos, name)
    195     eq_kwargs = {}
    197 if eq_file is not None:
--> 198     self.load_global_eq(eq_file, eq_type, **eq_kwargs)
    200 # Load global kinetics file if it exists
    201 if kinetics_kwargs is None:

File [~/programming/python/pyrokinetics/src/pyrokinetics/pyro.py:1504](http://localhost:8889/lab/tree/~/programming/python/pyrokinetics/src/pyrokinetics/pyro.py#line=1503), in Pyro.load_global_eq(self, eq_file, eq_type, **kwargs)
   1481 """
   1482 Reads a global equilibrium file, sets the property ``eq_file`` to that file
   1483 path, and sets the attribute ``eq`` to an Equilibrium.
   (...)
   1501     Equilibrium.
   1502 """
   1503 self.eq_file = eq_file  # property setter, converts to Path
-> 1504 self.eq = read_equilibrium(self.eq_file, eq_type, **kwargs)

File [~/programming/python/pyrokinetics/src/pyrokinetics/equilibrium/equilibrium.py:1054](http://localhost:8889/lab/tree/~/programming/python/pyrokinetics/src/pyrokinetics/equilibrium/equilibrium.py#line=1053), in read_equilibrium(path, file_type, **kwargs)
   1050 def read_equilibrium(
   1051     path: PathLike, file_type: Optional[str] = None, **kwargs
   1052 ) -> Equilibrium:
   1053     r"""A plain-function alternative to ``Equilibrium.from_file``."""
-> 1054     return Equilibrium.from_file(path, file_type=file_type, **kwargs)

File [~/programming/python/pyrokinetics/src/pyrokinetics/file_utils.py:246](http://localhost:8889/lab/tree/~/programming/python/pyrokinetics/src/pyrokinetics/file_utils.py#line=245), in ReadableFromFile.from_file(cls, path, file_type, **kwargs)
    244 # Infer reader type from path if not provided with file_type
    245 reader = cls._factory(str(path) if file_type is None else file_type)
--> 246 return reader(path, **kwargs)

File [~/programming/python/pyrokinetics/src/pyrokinetics/file_utils.py:93](http://localhost:8889/lab/tree/~/programming/python/pyrokinetics/src/pyrokinetics/file_utils.py#line=92), in AbstractFileReader.__call__(self, filename, *args, **kwargs)
     91 def __call__(self, filename: PathLike, *args, **kwargs) -> Any:
     92     """Forwards calls to ``read_from_file``."""
---> 93     return self.read_from_file(filename, *args, **kwargs)

File [~/programming/python/pyrokinetics/src/pyrokinetics/equilibrium/gacode.py:154](http://localhost:8889/lab/tree/~/programming/python/pyrokinetics/src/pyrokinetics/equilibrium/gacode.py#line=153), in EquilibriumReaderGACODE.read_from_file(self, filename, nR, nZ, clockwise_phi, cocos, neighbors)
    151 len_units = units.meter
    152 psi_units = units.weber [/](http://localhost:8889/) units.radian
--> 154 profiles = read_gacode_file(filename)
    156 psi = profiles.polflux * psi_units
    157 B_0 = profiles.bcentr * units.tesla

File [~/programming/python/pyrokinetics/src/pyrokinetics/equilibrium/gacode.py:63](http://localhost:8889/lab/tree/~/programming/python/pyrokinetics/src/pyrokinetics/equilibrium/gacode.py#line=62), in read_gacode_file(filename)
     61 # Check if relevant keys exist}
     62 if len(set(test_keys).intersection(data_dict.keys())) != len(test_keys):
---> 63     raise ValueError("EquilibriumReaderGACODE could not find all relevant keys")
     65 for key, value in data_dict.items():
     66     # If data has two columns, convert to 2D array
     67     setattr(data_object, key, np.squeeze(np.array(value)))

ValueError: EquilibriumReaderGACODE could not find all relevant keys

input.gacode.txt

[bpatel2107]

It looks like your input.gacode file is missing fpol $f(\psi)$, and this is one of the inputs that pyro checks for. We could probably get around this but ideally this data should be loaded. How did you generate the input.gacode file?

[jmcclena]

It came from running profiles_gen from a p-file and g-file. Weird, I see that fpol is in the OMFITtree, I guess OMFIT is deriving it from expro.

[jmcclena]

Hmmm... It actually does exist in the one generated from profiles_gen. OMFIT must drop it in the save or something.

[jmcclena]

@orso82 Do you remember why OMFIT is dropping fpol in the save? I assume is because its considered a derived quantity, but this quantity ends up noisier than the original profiles_gen value (and opposite sign?).

[bpatel2107]

@orso82 @jmcclena any updates on this?