diff --git a/docs/src/data/example-config.ts b/docs/src/data/example-config.ts index 38dc0ed..7e88c4e 100644 --- a/docs/src/data/example-config.ts +++ b/docs/src/data/example-config.ts @@ -212,8 +212,8 @@ export const exampleConfig: Record = { category: 'Stellarator particle orbits', setupTitle: 'Guiding centers in a generated GVEC stellarator', plotTitle: 'Guiding-center trajectories on 3D stellarator flux surfaces', plotAlt: 'Six guiding-center trajectories on nested magnetic flux surfaces of a five-field-period GVEC stellarator', figures: { 'flux-surfaces': { - alt: 'Nested GVEC flux surfaces in a poloidal cut beside pressure and rotational-transform profiles', - caption: 'A poloidal cut through the nested magnetic surfaces generated by GVEC. The pressure and rotational transform on the right are evaluated directly from the newly generated equilibrium state.', + alt: 'Rotational-transform and pressure profiles beside |B| on a poloidal cut through the nested GVEC flux surfaces', + caption: 'The rotational transform ι, with its lowest-order rational surface marked, and the pressure over the flux radius ρ, and |B| on the poloidal cut at ζ = 0 with the lines of constant ρ and θ: the nested magnetic surfaces generated by GVEC. All three are evaluated directly from the newly generated equilibrium state that Struphy uses.', }, 'poloidal-orbits': { alt: 'One panel per marker showing its guiding-center orbit projected on the poloidal plane', diff --git a/docs/src/examples/beltrami-sph.py b/docs/src/examples/beltrami-sph.py index 9bbf2b1..dc902eb 100644 --- a/docs/src/examples/beltrami-sph.py +++ b/docs/src/examples/beltrami-sph.py @@ -329,17 +329,8 @@ def error_trace(index, values, name, color): figure.update_xaxes(title_text="t", range=[0, float(times[-1])], row=1, col=2) figure.update_yaxes(title_text="relative error", type="log", range=[-7, -0.7], row=1, col=2) - last = len(times) - 1 - final_data = [ - contour, - marker_trace(last), - error_trace(last, velocity_error, "velocity RMS", "#00a884"), - error_trace(last, energy_error, "Hamiltonian drift", "#9b51e0"), - ] - # The frames name only the markers and errors (the contour stays), so build the final still. - still = go.Figure(data=final_data, layout=figure.layout) - still.layout.sliders[0].active = len(figure.frames) - 1 - save_figure(figure, "beltrami-sph", still=still, width=1100, height=680, show=show) + # The image shows the last frame, the end of the run (the contour stays as it is). + save_figure(figure, "beltrami-sph", frame=-1, width=1100, height=680, show=show) # The kernel reconstruction shows the simulated mass density independently of the marker view. # The exact divergence-free Beltrami transport preserves the initially uniform rho = 1. diff --git a/docs/src/examples/guiding-center-orbits.py b/docs/src/examples/guiding-center-orbits.py index 86261e8..980509d 100644 --- a/docs/src/examples/guiding-center-orbits.py +++ b/docs/src/examples/guiding-center-orbits.py @@ -283,11 +283,8 @@ def moving_traces(index): ], ) - # The still image shows the end of the run, with all trails. - # (A frame names only the moving traces, so the still is built from them and the fixed ones.) - still = go.Figure(data=list(figure.data[:first_moving]) + moving_traces(len(times) - 1), layout=figure.layout) - still.layout.sliders[0].active = len(frames) - 1 - save_figure(figure, "guiding-center-orbits", height=750, still=still, show=show) + # The image shows the last frame, the end of the run, with all trails. + save_figure(figure, "guiding-center-orbits", frame=-1, height=750, show=show) # One panel per particle: the orbit in the poloidal plane. panels = make_subplots(rows=2, cols=4, subplot_titles=labels, horizontal_spacing=0.03, vertical_spacing=0.12) diff --git a/docs/src/examples/gvec-equilibrium.py b/docs/src/examples/gvec-equilibrium.py index 015a0a5..4c005a0 100644 --- a/docs/src/examples/gvec-equilibrium.py +++ b/docs/src/examples/gvec-equilibrium.py @@ -40,6 +40,7 @@ maxwellians, ) from struphy.models import GuidingCenter +import plasma_plots from plasma_plots import save_figure # Keep both the GVEC solve and the following FEEC simulation deliberately @@ -371,85 +372,34 @@ def pproc(sim: Simulation, show: bool = False): ) save_figure(figure, "gvec-equilibrium", height=760, show=show) - # A poloidal slice reveals the nesting more quantitatively; the adjacent - # radial profiles come from the very same state file used by Struphy. - cut_radii = np.linspace(0.1, 1.0, 10) - # The cut at zeta = 0 for the flux-surface panel, and a few more across one field period: the + # The nesting of the flux surfaces, quantitatively: |B| on the poloidal plane at zeta = 0 with the + # flux-coordinate lines, beside the rotational transform (with its rational surfaces) and the + # pressure, all read from the very state file Struphy uses. + flux = plasma_plots.from_gvec( + equilibrium.state.evaluate("mod_B", "pos", "iota", "p", "theta_P", "N_FP", rho=17, theta=64, zeta=40) + ) + # The cut comes last, so that its colour bar sits at the right edge of the figure. + with plasma_plots.figure(1, 3, backend="plotly", title="GVEC equilibrium profiles and flux geometry") as profiles: + flux.iota.plasma.plot.lineout(rationals=4, title="rotational transform ι", ax=profiles[0]) + flux.p.plasma.plot.lineout(title="pressure p", ax=profiles[1]) + flux.mod_B.plasma.plot.slice( + coords="physical", + plane="RZ", + zeta=0.0, + overlays={"coordinate_lines": {"rho": 5, "theta_P": 8}}, + title="|B| and flux surfaces (ζ = 0)", + ax=profiles[2], + ) + save_figure(profiles, "gvec-equilibrium-flux-surfaces", width=1300, show=show) + + # The cut at zeta = 0 and a few more across one field period, for the orbit panels: the # cross-section of a stellarator turns with the toroidal angle, so a projected orbit lives inside # the envelope of all of them rather than on any single cut. + cut_radii = np.linspace(0.1, 1.0, 10) panel_angles = np.linspace(0.0, 2.0 * np.pi / int(equilibrium.state.nfp), 5) - cut_data = equilibrium.state.evaluate( - "pos", rho=cut_radii, theta=181, zeta=np.append(np.array((0.0,)), panel_angles) - ) - cut_positions = np.asarray(cut_data["pos"])[..., :1] - panel_positions = np.asarray(cut_data["pos"])[..., 1:] - profile_radii = np.linspace(0.0, 1.0, 101) - profile_data = equilibrium.state.evaluate( - "iota", "p", rho=profile_radii, theta=0, zeta=0 - ) - - profiles = make_subplots( - rows=1, - cols=2, - specs=[[{}, {"secondary_y": True}]], - subplot_titles=( - "Poloidal flux surfaces (ζ = 0)", - "Radial equilibrium profiles", - ), - horizontal_spacing=0.14, - ) - for index, radius in enumerate(cut_radii): - # Not x, y, z: those hold the orbits, which the panels below still need. - cut_x = cut_positions[0, index, :, 0] - cut_y = cut_positions[1, index, :, 0] - cut_z = cut_positions[2, index, :, 0] - radial_position = np.sqrt(cut_x**2 + cut_y**2) - profiles.add_scatter( - x=np.append(radial_position, radial_position[0]), - y=np.append(cut_z, cut_z[0]), - mode="lines", - line={"color": "#168aad", "width": 1.5 + 1.2 * radius}, - opacity=0.35 + 0.65 * radius, - name=f"ρ = {radius:.1f}", - legendgroup="surfaces", - showlegend=index in (0, len(cut_radii) - 1), - row=1, - col=1, - ) - profiles.add_scatter( - x=profile_radii, - y=np.asarray(profile_data["p"]), - mode="lines", - name="pressure p", - line={"color": "#d62828", "width": 3}, - row=1, - col=2, - secondary_y=False, - ) - profiles.add_scatter( - x=profile_radii, - y=np.asarray(profile_data["iota"]), - mode="lines", - name="rotational transform ι", - line={"color": "#f4a261", "width": 3}, - row=1, - col=2, - secondary_y=True, - ) - profiles.update_xaxes(title_text="R", scaleanchor="y", scaleratio=1, row=1, col=1) - profiles.update_yaxes(title_text="Z", row=1, col=1) - profiles.update_xaxes(title_text="normalized flux radius ρ", row=1, col=2) - profiles.update_yaxes(title_text="pressure p", row=1, col=2, secondary_y=False) - profiles.update_yaxes( - title_text="rotational transform ι", row=1, col=2, secondary_y=True - ) - profiles.update_layout( - title="GVEC flux geometry and equilibrium profiles", - template="plotly_white", - legend={"orientation": "h", "y": -0.2}, - margin={"l": 70, "r": 70, "t": 90, "b": 100}, + panel_positions = np.asarray( + equilibrium.state.evaluate("pos", rho=cut_radii, theta=181, zeta=panel_angles)["pos"] ) - save_figure(profiles, "gvec-equilibrium-flux-surfaces", show=show) # One panel per marker, the orbit projected on the (R, Z) plane, as in the tokamak example. The # cross-section of a stellarator turns with the toroidal angle, so the surfaces drawn behind each diff --git a/submodules/plasma-plots b/submodules/plasma-plots index 9bbc3c1..43dbd03 160000 --- a/submodules/plasma-plots +++ b/submodules/plasma-plots @@ -1 +1 @@ -Subproject commit 9bbc3c175f17037517de012b6ccad5031dbd6477 +Subproject commit 43dbd032cf53be6649e87eda5a66478f197e045c