Analysis

%load_ext autoreload
%autoreload 2

import matplotlib.pyplot as plt
import yt
from rich import print
from icecream import ic
import panel as pn
from yt.data_objects.time_series import SimulationTimeSeries
from yt.data_objects.static_output import Dataset

from utils.plot import *
from main import setup

dim = 1
beta = 0.25
theta = 60.0
eta = 10.0

meta = setup(dim, beta, theta, eta)

Fields

from utils import load_ts_all
ts_field, ts_part = load_ts_all(meta) 

ic(len(ts_field))
ds_field: Dataset = ts_field[-1]
ds_part: Dataset = ts_part[0]

ic| len(ts_field): 101

ds_part.field_list

[('all', 'particle_cpu'),
 ('all', 'particle_id'),
 ('all', 'particle_momentum_x'),
 ('all', 'particle_momentum_y'),
 ('all', 'particle_momentum_z'),
 ('all', 'particle_position_x'),
 ('all', 'particle_weight'),
 ('ions', 'particle_cpu'),
 ('ions', 'particle_id'),
 ('ions', 'particle_momentum_x'),
 ('ions', 'particle_momentum_y'),
 ('ions', 'particle_momentum_z'),
 ('ions', 'particle_position_x'),
 ('ions', 'particle_weight'),
 ('nbody', 'particle_cpu'),
 ('nbody', 'particle_id'),
 ('nbody', 'particle_momentum_x'),
 ('nbody', 'particle_momentum_y'),
 ('nbody', 'particle_momentum_z'),
 ('nbody', 'particle_position_x'),
 ('nbody', 'particle_weight')]

ps = "ions"
_part_weight_field = (ps, "particle_weight")

p = yt.ParticlePlot(
    ds_part, ('ions', 'particle_position_x'), [('ions', 'particle_momentum_y'), ("ions", "particle_momentum_z")],
    z_fields=_part_weight_field
)
p.set_log(field = 'all', log=False)

Archive

yt

ts: SimulationTimeSeries = yt.load('diags/diag1??????')
# ts = yt.load('./diags/diag???0032')

def plot(ds, normalize = True, **kwargs):
    ad = ds.all_data()
    fields = ["Bx", "By", "Bz"]

    match meta["dim"]:
        case 1: pos = "x"
        case 2: pos = "y"
        
    pos = ad[pos]
    
    if normalize:
        pos = pos / meta['d_i']
    
    for field in fields:
        plt.plot(pos, ad[field], label=field, **kwargs)
    
    plt.xlabel("x ($d_i$)")

for i, part_df in enumerate(ts):
    alpha = (i + 1) / (len(ts)+1)
    plot(part_df, alpha=alpha)
    plt.show()  # Show each plot separately

i = 4
_ts = ts[0:i]
for i, part_df in enumerate(_ts):
    alpha = (i + 1) / (len(_ts)+1)
    plot(part_df, alpha=alpha)
    plt.show()  # Show each plot separately

yt.SlicePlot(part_df, "z", ("boxlib", "Bz"))

part_df.all_data()

fields = [
    ("Bx"),
    ("By"),
    ("Bz"),
    ("mesh", "magnetic_field_strength"),
]

for part_df in ts.piter():
    p = yt.plot_2d(part_df, fields=fields)
    p.set_log(fields, False)
    fig = p.export_to_mpl_figure((2, 2))
    fig.tight_layout()
    fig.savefig(f"figures/{part_df}_magnetic_field.png")

Average magnetic field

def plot_avg(ds):
    fields = [
        ("Bx"),
        ("By"),
        ("Bz"),
    ]

    ad = ds.all_data()
    df = ad.to_dataframe(fields + ["x", zaxis])
    # compute the magnetic field strength
    df = df.assign(B=lambda x: (x.Bx**2 + x.By**2 + x.Bz**2) ** 0.5)

    axes = df.groupby(zaxis).mean().plot(y=fields + ["B"], subplots=True)
    return axes[0].figure

def plot_avg_ts(i):
    return plot_avg(ts[i])
    
time_widget = pn.widgets.IntSlider(name="Time", value=1, start=0, end=len(ts)-1)
bound_plot = pn.bind(plot_avg_ts, i=time_widget)

pn.Column(time_widget, bound_plot)

NameError: name 'ts' is not defined

for part_df in ts.piter():
    plot_avg(part_df)

part_df.print_stats()
print(part_df.field_list)

grid = part_df.r[:,:,:]
obj = grid.to_xarray(fields=fields)