PlasmaControl / DESC

Stellarator Equilibrium and Optimization Suite
MIT License
87 stars 22 forks source link

TypeError in 'plot_comparison' #1202

Closed Taegeun-Jeong closed 2 weeks ago

Taegeun-Jeong commented 3 weeks ago

Hello,

When I run the plot_comparison function, I get some errors. It seems that it is necessary for the user to define "phi" parameters to avoid this TypeError.

TypeError: reshape total size must be unchanged, got new_sizes (5, 1, 180) for shape (1080,).

import numpy as np
import desc.io

from desc.geometry import FourierRZToroidalSurface
from desc.plotting import plot_comparison

# Load Equilibrium and Make Surface

eq0 = desc.io.load("QA.h5")
eq1 = desc.io.load("QA3_1_iota_0_42_AR5.h5")

surf = FourierRZToroidalSurface(
    R_lmn =   [0.97, 0.33],
    modes_R = [(0, 0), (-1, 0)],
    Z_lmn =   [0.4],
    modes_Z = [(1, 0)],
)

# No error occurs for some Equilibria and Surfaces
plot_comparison([eq0, eq1], rho=np.array(1.0), theta=0);
plot_comparison([surf, eq0], rho=np.array(1.0), theta=0);

# However, error occurs for some Equilibrium and Surface.
plot_comparison([surf, eq1], rho=np.array(1.0), theta=0);

# This error no longer occurs, when I specify the phi parameter.
plot_comparison([surf, eq1], rho=np.array(1.0), theta=0, phi=np.linspace(0, 2*np.pi/eq0.NFP * 5/6, 6));

Attached is the Jupyter notebook and *.h5 file I used. Bug_report_plot_comparison.zip

YigitElma commented 3 weeks ago

So, I guess the main use case of plot_comparison function is for equilibria with the same number of field periods. However, you have a surface with NFP=1, eq0.NFP=2 and eq1.NFP=3. So, although you don't get an error for 2 equilibria case, in fact what you plot is wrong (this is the case if you want to use NFP of each equilibrium as the titles suggest, but if you want to plot the surfaces at same phi angle without even considering NFP) because inside plot_comparison function (if you don't give phi) NFP is set using the first objects NFP. For example, here is what I think you want to plot,

fig, ax = plot_comparison([eq0], rho=np.array(1.0), theta=0, labels=["eq0"], color="red");
plot_comparison([eq1], rho=np.array(1.0), theta=0, labels=["eq1"], ax=ax);

image

We should throw an error for inputs with different NFP.

YigitElma commented 3 weeks ago

For the surface and eq1 case, if you want to able to draw that correctly, you can use this,

plot_comparison([surf, eq1], rho=1., theta=0, phi=np.linspace(0, 2*np.pi/eq1.NFP, 6, endpoint=False));

image But this wouldn't be correct if your surface wasn't axisymmetric!

dpanici commented 3 weeks ago

1204 should explain why these errors occur (in one case they are a bug, in the other it is an unforeseen use case that is ambiguous) and resolve them.

If you still want to plot multiple equilibria which have differing NFPs where they are both nonaxisymmetric, you should instead make multiple calls to plot_surfaces and pass in the desired phi values to ensure that the resulting plot is comparing what you expect, like

phi = np.linspace(0,2*np.pi,6,endpoint=False)
fig,ax = plot_surfaces(eq0, phi=phi)
plot_surfaces(eq1, phi=phi,ax=ax)