dstansby / pfsspy

Potential Field Source Surface model package for Python
https://pfsspy.readthedocs.io/
Other
40 stars 17 forks source link

Exception when tracing fieldlines in pfsspy v1.0.0 #330

Closed wtbarnes closed 2 years ago

wtbarnes commented 2 years ago

When running the fieldline tracing in this notebook https://github.com/wtbarnes/hinode-iris-2021-sunpy-tutorial, I get the following exception:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
/var/folders/l0/9klzdmq17mxdq3r_00f903s0000542/T/ipykernel_32715/397003646.py in <module>
----> 1 fieldlines = tracer.trace(SkyCoord(seeds_eis), output,)

~/miniconda3/envs/sunpy-hinode-iris-2021/lib/python3.9/site-packages/pfsspy/tracing.py in trace(self, seeds, output)
    182 
    183         xs = [np.stack(pfsspy.coords.strum2cart(x[:, 2], x[:, 1], x[:, 0]), axis=-1) for x in xs]
--> 184         flines = [fieldline.FieldLine(x[:, 0], x[:, 1], x[:, 2], output) for x in xs]
    185         return fieldline.FieldLines(flines)
    186 

~/miniconda3/envs/sunpy-hinode-iris-2021/lib/python3.9/site-packages/pfsspy/tracing.py in <listcomp>(.0)
    182 
    183         xs = [np.stack(pfsspy.coords.strum2cart(x[:, 2], x[:, 1], x[:, 0]), axis=-1) for x in xs]
--> 184         flines = [fieldline.FieldLine(x[:, 0], x[:, 1], x[:, 2], output) for x in xs]
    185         return fieldline.FieldLines(flines)
    186 

~/miniconda3/envs/sunpy-hinode-iris-2021/lib/python3.9/site-packages/pfsspy/fieldline.py in __init__(self, x, y, z, output)
    131         # the source surface
    132         atol = 0.1
--> 133         self._is_open = np.abs(self._r[0] - self._r[-1]) > atol
    134         self._polarity = -np.sign(self._r[0] - self._r[-1]) * self._is_open
    135 

IndexError: index 0 is out of bounds for axis 0 with size 0

Notably, this does not seem to happen with version 0.6.6.

I will try to follow up with a much more pared down MWE. I realized that one is overly complicated. I just wanted to preserve this exception as a reminder.

dstansby commented 2 years ago

Are you using a custom step size in the tracer? If so, what is it?

wtbarnes commented 2 years ago

I am not. This is just using the FORTRAN tracer with no additional Kwargs

STBadman commented 2 years ago

I am also seeing this issue in 1.0.0 - I get it when downloading and running the example as well as in my own code,

Strangely, i don't get the error if i restrict it to a single set of longitudes [0,360] at 0 latitude and start the tracing from the source surface, but every other case I've tried gives the same error.

dstansby commented 2 years ago

Thanks for providing a (link to) a self contained example - I'll take a look at this next week. An initial guess is it's issues with field lines traced from near the poles, which would explain why a set of seeds at the equator is fine.

wtbarnes commented 2 years ago

FWIW, in the example application I'm seeing this exception in, I'm limiting my seed points to a single AR that is very far from the poles.

dstansby commented 2 years ago

What version of streamtracer are you (both) using? Using streamtracer 1.1.1 the open/closed map example seems to work fine for me.

STBadman commented 2 years ago

I have:

Python 3.8 pfsspy 1.0.0 Streamtracer 1.1.1 Sunpy 3.0.2 astropy 4.2.1

And get the error by downloading the example .py file locally and running >>> python open_closed_map.py

dstansby commented 2 years ago

Thanks - I think I've found the issue and the solution.