PlasmaControl / DESC

Stellarator Equilibrium and Optimization Suite
MIT License
96 stars 26 forks source link

Advice on plotting #1343

Open mohawk811 opened 1 day ago

mohawk811 commented 1 day ago

Hello!

I was wondering if I could get some input on some custom plots I'm trying to do. For this code I want to plot QS error on the Y axis and aspect ratio on the X axis. I expect to see that as I go to higher aspect ratio QS error should decrease but my plot is showing otherwise. Could I get some advice or input on how to show this?

import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 22})
import numpy as np

# Lists to store aspect ratios, elongations, QS errors, and corresponding 'a' values
aspect_ratios = []
elongations = []
qs_errors = []
file_aspect = []

# Loop through each equilibrium
for a in range(4, 20):
    # Generate the filename based on the current elongation value 'a'
    filename = f'eq_opt_qs_3_aspect_{a}_10_29.h5'

    # Load the equilibrium object
    eq = desc.io.load(filename)

    # Define the grid
    grid = LinearGrid(M=eq.M_grid, N=eq.N_grid, NFP=eq.NFP, rho=np.array([.2,.6,1.0]), sym=True)

    # Compute the aspect ratio R0/a
    aspect_ratio = eq.compute('R0/a', grid=grid)['R0/a'].max()

    # Compute the elongation a_major/a_minor
    elongation_value = eq.compute('a_major/a_minor', grid=grid)['a_major/a_minor'].max()

    # Compute the QS error
    qs_error = eq.compute('f_C', grid=grid, helicity=(0,1))['f_C'].max()  # Assuming 'f_C' is the QS error

    # Append the values to the lists
    aspect_ratios.append(aspect_ratio)
    elongations.append(elongation_value)
    qs_errors.append(qs_error)
    file_aspect.append(a)  # Store the value of 'a' from the filename

# Plotting Aspect Ratio vs. Elongation with QS Error Color Scale
plt.figure(figsize=(10, 6))
sc = plt.scatter(aspect_ratios, elongations, c=qs_errors, cmap='viridis', marker='o', vmin=0, vmax=.5)

for i in range(len(aspect_ratios)):
    # Slight offset for the text position
    x_offset = 0.02 * (i % 2 - 0.5)
    y_offset = 0.02 * (i % 2 - 0.5)

    # Aspect ratio and elongation annotation
    plt.text(aspect_ratios[i] + x_offset, elongations[i] + y_offset,
             f'({aspect_ratios[i]:.2f}, {elongations[i]:.2f})',
             fontsize=9, ha='right', va='bottom', bbox=dict(facecolor='white', alpha=0.5))

    # Annotate with the 'a' value from the filename in red with a slight offset
    plt.text(aspect_ratios[i] - x_offset, elongations[i] - y_offset,
             f'a={file_aspect[i]:.2f}', fontsize=9, ha='left', va='top', color='red',
             bbox=dict(facecolor='white', alpha=0.5))

plt.xlabel('Aspect Ratio (R0/a)')
plt.ylabel('Average Elongation (a_major/a_minor)')
plt.title('Aspect Ratio vs. Average Elongation with QS Error Color Scale')
plt.colorbar(sc, label='QS Error (f_C)')  # Add a colorbar to represent QS error
plt.grid(True)
plt.show()

# New plot: Aspect Ratio vs. QS Error (f_C)
plt.figure(figsize=(10, 6))
plt.plot(aspect_ratios, qs_errors, marker='o', color='b', label='QS Error (f_C)')

# Annotating 'a' values
for i in range(len(aspect_ratios)):
    plt.text(aspect_ratios[i], qs_errors[i], f'a={file_aspect[i]}', fontsize=9, ha='right', va='bottom', color='red')

plt.xlabel('Aspect Ratio (R0/a)')
plt.ylabel('QS Error (f_C)')
plt.title('Aspect Ratio vs. QS Error')
plt.grid(True)
plt.legend()
plt.show()
mohawk811 commented 1 day ago

Here are the generated plots:

Screenshot 2024-10-31 at 1 33 02 PM
mohawk811 commented 1 day ago
Screenshot 2024-10-31 at 1 32 53 PM
dpanici commented 16 hours ago

hard to say much without knowing the specifics of the optimization and how good of QP solutions the individual equilibria are, but some ideas for the plots