freegs-plasma / freegs

Free boundary Grad-Shafranov solver
http://freegs.readthedocs.io/en/latest/
GNU Lesser General Public License v3.0
91 stars 54 forks source link

Secondary x-points discarded by lack of flux-monotonicity #55

Closed pabloprf closed 2 years ago

pabloprf commented 3 years ago

Starting in line 196 in critical.py, x-points that feature a non-monotonic flux in the line that connects them to the o-points are discarded. Problem is that secondary x-points (for example, for x-point target divertor) are discarded. Like the one in the figure attached.

Screen Shot 2021-03-03 at 10 54 28 AM

My questions are:

  1. Why is the discard_xpoints boolean commented in that routine?
  2. What is the goal of this discarding function? (sorry if this is too naive)
  3. Is there a suggested workaround for secondary x-points?

Thanks!

bendudson commented 3 years ago

Hi @pabloprf It's been a while so I may be forgetting something, but I think in this order:

  1. The goal is to identify the primary X-point, so that the normalised psi can be calculated correctly. The rest of the code (at least as it is in Github) doesn't care about the other X-points found by this routine (the X-point control system constraints are independent of this)

  2. Not sure why; this may have been a temporary thing that got accidentally committed

  3. As long as the routine can identify the primary X-point, then the non-monotonic X-points don't need to be discarded.

Perhaps these X-points could be put into a separate list here: https://github.com/bendudson/freegs/blob/master/freegs/critical.py#L226

so that after the X-points are sorted here https://github.com/bendudson/freegs/blob/master/freegs/critical.py#L232 this non-monotonic list can be appended to the xpoint list. That way they would come after the primary X-point (which is the first in xpoint).

pabloprf commented 3 years ago

Thanks @bendudson. It is useful to know that only the primary x-point is used in the code. Actually one quick fix to this would be to uncomment discard_xpoints (and have the keyword argument discard_xpoints=True). This way, the iteration routines are correctly using the primary x-point. Then, when the user wants to get all the x-points (e.g. for plotting at the end), we simply use find_critical(R, Z, psi, discard_xpoints=False).

pabloprf commented 3 years ago

A related question is that if I run the entire workflow without the discarding function (keeping all x-points), I get very different final equilibria. As we have discussed in the past, I'm running FreeGS without control constrains (coil currents completely fixed, not allowed to vary). Do you know why that may be?

bendudson commented 3 years ago

Unfortunately if you don't discard those equilibria, then when the X-points are sorted by psi (https://github.com/bendudson/freegs/blob/master/freegs/critical.py#L232) it will pick the wrong X-point as the primary. That will lead the poloidal flux normalisation to be different, change current profiles etc. That's why I think you have to separate the non-monotonic ones, sort the ones which are monotonic, and then put the non-monotonic ones at the end. That way they're not included in the sorting.

pabloprf commented 3 years ago

I'm finding an interesting behavior. Without discarding, even if I force the first x-point to be the primary (for example, by sorting by R, and I'm making sure the primary is always first), the existence of the non-monotonic x-points at the end of the x-points list causes the solver to find very different equilibria. The only thing that comes to mind that could be causing this is that the x-points are again sorted by psi somewhere outside of find_critical.

bendudson commented 3 years ago

That is curious. I suspect you must be right, that they're being sorted somewhere else.

pabloprf commented 3 years ago

Just FYI. I found out that those non-monotonic x-points are used here: https://github.com/bendudson/freegs/blob/9dfc4b7201263fee3385b9bfc67aae167564c9a8/freegs/critical.py#L276

If I change that line to something like for rx, zx, _ in xpoint[:2]:, then the problem disappears and the final equilibria is the same as one where the non-monotonic x-points have been discarded. In any case, I will run it as I said here: https://github.com/bendudson/freegs/issues/55#issuecomment-789865118

bendudson commented 3 years ago

Ah good find! Difficult... That loop needs to go over the primary X-point, and any others which are sufficiently close. In a connected double null plasma, both upper and lower X-points, but in a single null just one.

pabloprf commented 3 years ago

Sounds good, thanks @bendudson. Do you want me to submit a PR uncommenting discard_xpoints and changing the default to discard_xpoints=True? The exact same behavior we have until now is recovered, but now the user could also retrieve all x-points if needed.

bendudson commented 3 years ago

Thanks @pabloprf that would be great when you have time. I think I see how to implement a fix, but I don't have a good test case. It would be nice to have a unit test to check the result of this function with secondary X-points, but I don't see me having time to do it for a few weeks.

pabloprf commented 3 years ago

Thanks @bendudson. I have created a PR to add the temporary fix. If you're interested I can put together a test case so that you can try your other fix when you have time.