Post-processing and standard plots#
This tutorial introduces the standardized post-processing interface. We run a small Vlasov–Ampère example, get its output as an autocomplete-friendly Output, and make the plots most commonly used to inspect a simulation.
For a production run you can skip the simulation setup and open its output folder with struphy.Output("path/to/sim") instead.
[1]:
import os
import tempfile
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import HTML
from struphy import (
BinningPlot,
BoundaryParameters,
ButcherTableau,
DerhamOptions,
EnvironmentOptions,
KernelDensityPlot,
LoadingParameters,
SavingParameters,
Simulation,
SortingParameters,
Time,
WeightsParameters,
domains,
equils,
grids,
maxwellians,
perturbations,
)
from struphy.models import Maxwell, ViscousEulerSPH, VlasovAmpereOneSpecies
Create a compact demonstration run#
Post-processing operates on a completed run. The small setup below saves an electric field, a few marker trajectories, scalar diagnostics, and a binned \((\eta_1,v_1)\) distribution. These are the main output types handled by the plotting interface.
[2]:
def build_model():
model = VlasovAmpereOneSpecies(alpha=1.0, epsilon=-1.0, with_B0=False)
model.em_fields.e_field.save_data = True
model.em_fields.phi.save_data = True
model.kinetic_ions.var.save_data = True
model.propagators.push_eta.options = model.propagators.push_eta.Options()
model.propagators.coupling_va.options = model.propagators.coupling_va.Options()
model.initial_poisson.options = model.initial_poisson.Options(stab_mat="M0")
binplot = BinningPlot(
slice="e1_v1",
n_bins=(32, 32),
ranges=((0.0, 1.0), (-5.0, 5.0)),
)
model.kinetic_ions.set_markers(
loading_params=LoadingParameters(ppc=32, seed=1234),
weights_params=WeightsParameters(control_variate=True),
boundary_params=BoundaryParameters(),
sorting_params=SortingParameters(boxes_per_dim=(4, 1, 1), do_sort=True),
saving_params=SavingParameters(n_markers=12, binning_plots=(binplot,)),
)
background = maxwellians.Maxwellian3D(n=(1.0, None))
model.kinetic_ions.var.add_background(background)
density_mode = perturbations.ModesCos(ls=(1,), amps=(1e-3,))
model.kinetic_ions.var.add_initial_condition(maxwellians.Maxwellian3D(n=(1.0, density_mode)))
return model
model = build_model()
/opt/hostedtoolcache/Python/3.10.21/x64/lib/python3.10/site-packages/struphy/models/species.py:215: UserWarning: Override equation parameter self.alpha =1.0
warnings.warn(f"Override equation parameter {self.alpha =}")
/opt/hostedtoolcache/Python/3.10.21/x64/lib/python3.10/site-packages/struphy/models/species.py:222: UserWarning: Override equation parameter self.epsilon =-1.0
warnings.warn(f"Override equation parameter {self.epsilon =}")
[3]:
demo_tmp = tempfile.TemporaryDirectory(prefix="struphy_postprocessing_")
demo_root = demo_tmp.name
env = EnvironmentOptions(
out_folders=demo_root,
sim_folder="vlasov_ampere_demo",
save_restart=False,
)
sim = Simulation(
model=model,
env=env,
time_opts=Time(dt=0.05, Tend=5.0),
domain=domains.Cuboid(r1=2 * 3.141592653589793),
equil=equils.HomogenSlab(),
grid=grids.TensorProductGrid(num_elements=(16, 1, 1)),
derham_opts=DerhamOptions(degree=(2, 1, 1)),
)
out = sim.run(profiling_activated=True)
print(f"Raw output: {sim.env.path_out}")
Stabilizing Poisson solve with self.options.sigma_1 =1e-14
Time stepping: 100%|██████████| 100/100 [00:01<00:00, 86.06step/s]
╭────────────────────────────────────────────────────────────────────────────────────────────╮
│ region % session total [s] │
├────────────────────────────────────────────────────────────────────────────────────────────┤
│ scope_profiler.session 100.00% 1.319101 │
│ └─ (own) 1.15% 0.015221 │
│ └─ setup: total 11.69% 0.154184 │
│ │ └─ (own) 0.05% 0.000686 │
│ │ └─ setup: allocate 5.80% 0.076512 │
│ │ │ └─ (own) 0.01% 0.000070 │
│ │ │ └─ setup: feec 3.42% 0.045133 │
│ │ │ │ └─ (own) 0.01% 0.000084 │
│ │ │ │ └─ setup: derham 3.41% 0.044998 │
│ │ │ │ └─ setup: mass ops 0.00% 0.000009 │
│ │ │ │ └─ setup: basis ops 0.00% 0.000008 │
│ │ │ │ └─ setup: projected equil 0.00% 0.000034 │
│ │ │ └─ setup: variables 0.24% 0.003157 │
│ │ │ │ └─ (own) 0.01% 0.000139 │
│ │ │ │ └─ setup var: em_fields.e_field 0.01% 0.000181 │
│ │ │ │ └─ setup var: em_fields.phi 0.00% 0.000046 │
│ │ │ │ └─ setup var: kinetic_ions.var 0.21% 0.002792 │
│ │ │ │ │ └─ (own) 0.20% 0.002602 │
│ │ │ │ │ └─ do_sort 0.01% 0.000190 │
│ │ │ │ │ │ └─ (own) 0.00% 0.000042 │
│ │ │ │ │ │ └─ put_particles_in_boxes 0.01% 0.000147 │
│ │ │ └─ setup: propagators 1.51% 0.019877 │
│ │ │ │ └─ (own) 0.00% 0.000046 │
│ │ │ │ └─ setup prop: PushEta 0.00% 0.000040 │
│ │ │ │ └─ setup prop: VlasovAmpereCoupling 1.50% 0.019791 │
│ │ │ └─ setup: helpers 0.63% 0.008276 │
│ │ │ │ └─ (own) 0.28% 0.003674 │
│ │ │ │ └─ accum: charge_density_0form 0.02% 0.000228 │
│ │ │ │ │ └─ (own) 0.00% 0.000029 │
│ │ │ │ │ └─ kernel: charge_density_0form 0.01% 0.000142 │
│ │ │ │ │ └─ accum comm: charge_density_0form 0.00% 0.000057 │
│ │ │ │ └─ solve: PoissonSolve 0.33% 0.004319 │
│ │ │ │ └─ update_feec_variables 0.00% 0.000055 │
│ │ └─ setup: run metadata 0.08% 0.001120 │
│ │ └─ setup: data storage 0.42% 0.005550 │
│ │ └─ setup: geometry vtk 4.28% 0.056405 │
│ │ └─ setup: plasma params 0.23% 0.002981 │
│ │ └─ setup: initial diagnostics 0.12% 0.001599 │
│ │ └─ setup: hdf5 datasets 0.71% 0.009329 │
│ └─ model.integrate (100x) 41.58% 0.548451 │
│ │ └─ (own) 0.08% 0.001086 │
│ │ └─ prop: PushEta (100x) 4.96% 0.065436 │
│ │ │ └─ (own) 2.40% 0.031652 │
│ │ │ └─ pusher: push_eta_stage (100x) 2.56% 0.033783 │
│ │ │ │ └─ (own) 1.12% 0.014721 │
│ │ │ │ └─ kernel: push_eta_stage (400x) 1.45% 0.019063 │
│ │ └─ prop: VlasovAmpereCoupling (100x) 36.53% 0.481930 │
│ │ │ └─ (own) 3.17% 0.041809 │
│ │ │ └─ accum: vlasov_maxwell (100x) 10.95% 0.144435 │
│ │ │ │ └─ (own) 0.49% 0.006505 │
│ │ │ │ └─ kernel: vlasov_maxwell (100x) 5.13% 0.067735 │
│ │ │ │ └─ accum comm: vlasov_maxwell (200x) 5.32% 0.070195 │
│ │ │ └─ solve: SchurSolver (100x) 19.54% 0.257770 │
│ │ │ └─ pusher: push_v_with_efield (100x) 2.01% 0.026484 │
│ │ │ │ └─ (own) 0.61% 0.008086 │
│ │ │ │ └─ kernel: push_v_with_efield (100x) 1.39% 0.018398 │
│ │ │ └─ update_feec_variables (100x) 0.87% 0.011431 │
│ └─ diagnostics (100x) 5.47% 0.072151 │
│ └─ save data (101x) 40.11% 0.529094 │
╰────────────────────────────────────────────────────────────────────────────────────────────╯
╭─ Info ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────╮
│ Summary: ../../../../../../../../tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_demo/profiling_data.h5 (1 rank) │
│ │
│ Explore: │
│ Inspect: scope-profiler inspect ../../../../../../../../tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_demo/profiling_data.h5 │
│ TUI: scope-profiler tui ../../../../../../../../tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_demo/profiling_data.h5 │
│ │
│ Visualize and export: │
│ Plot: scope-profiler plot default ../../../../../../../../tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_demo/profiling_data.h5 -o plots --show │
│ Report: scope-profiler report ../../../../../../../../tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_demo/profiling_data.h5 -o report.html │
│ Export: scope-profiler export plot-data ../../../../../../../../tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_demo/profiling_data.h5 -o data │
│ Lines: scope-profiler line-profile ../../../../../../../../tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_demo/profiling_data.h5 │
│ │
│ Compare runs: │
│ Diff: scope-profiler diff BASE.h5 CANDIDATE.h5 │
│ Check: scope-profiler check BASE.h5 CANDIDATE.h5 │
│ │
│ Durations are in seconds. │
│ Regions may nest, so the summed total can exceed the wall-clock time. │
│ % session uses wall-clock coverage; overlapping recursive calls count once. │
╰─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────╯
Raw output: /tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_demo
Process and load the output#
sim.run() returns an Output, which is also available later as sim.output. Scalars are read directly from the raw output; fields and particle products need post-processing, which runs with default options the first time they are accessed.
To choose options, call out.pproc() first. Output evaluates saved FEEC fields and organizes particle diagnostics; physical=True additionally creates physical field components. Existing products made with the same options are reused, so re-running a cell is cheap.
Individual products are standard xarray.DataArray objects with named dimensions, coordinates, units, and labels. Time is in Struphy units, in which the models’ analytic results are written; seconds come along as the coordinate t_seconds, and out.with_time_units("physical") gives an independent view with t in seconds. Arrays are loaded only when accessed. The saved configuration is available through out.domain, out.model, and out.time_opts; no simulation object is
created.
[4]:
out.pproc(physical=True)
Post-processing path /tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_demo
Reading hdf5 data of following species:
em_fields:
e_field: <HDF5 group "/feec/em_fields/e_field" (3 members)>
phi: <HDF5 dataset "phi": shape (101, 20, 3, 3), type "<f8">
Creation of Struphy Fields done.
Evaluating fields ...
100%|██████████| 101/101 [00:00<00:00, 482.04it/s]
Evaluation of 12 marker orbits for kinetic_ions
100%|██████████| 101/101 [00:00<00:00, 848.04it/s]
Evaluation of distribution functions for kinetic_ions
0 starting post-processing of distribution functions for /tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_demo/post_processing/kinetic_data/kinetic_ions ...
100%|██████████| 1/1 [00:00<00:00, 1933.75it/s]
0%| | 0/1 [00:00<?, ?it/s]rank = 0 ----------------------------
self._pproc_rank =0 with xp.sum(data) =np.float64(10342.800145647721) and xp.sum(data_df) =np.float64(102.7174290693037)
self._pproc_rank =0 with xp.sum(data) =np.float64(10342.800145647721) and xp.sum(data_df) =np.float64(102.7174290693037)
self._pproc_rank =0 done.
100%|██████████| 1/1 [00:00<00:00, 29.08it/s]
[4]:
Output('/tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_demo', processed=True)
Products sit on the run under the species that produced them, so VS Code and interactive shells complete them as you type: out.kinetic_ions.e1_v1_density.f, out.kinetic_ions.orbits, out.em_fields.phi_log. To discover the products a particular namespace contains, inspect its flat, lazy .catalog: list(out.kinetic_ions.catalog). The grouped views out.fields, out.distributions, out.densities and out.orbits show the same products by kind. In scripts,
out.evaluate("kinetic_ions/f", dataset="e1_v1_density/f") looks up one explicitly.
[5]:
out.info()
Output: /tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_demo
Key Description Load with
---------------------------------- ------------------------------------------------- ---------
electric_energy scalar time series (electric energy) out.evaluate('scalars', variables='electric_energy')
em_fields/e_field field (e_field) out.evaluate('em_fields/e_field')
em_fields/e_field_xyz field (e_field_xyz) out.evaluate('em_fields/e_field_xyz')
em_fields/phi field (phi) out.evaluate('em_fields/phi')
em_fields/phi_xyz field (phi_xyz) out.evaluate('em_fields/phi_xyz')
kinetic_energy scalar time series (kinetic energy) out.evaluate('scalars', variables='kinetic_energy')
kinetic_ions marker trajectories (x, y, z, v1, v2, v3, weight) out.evaluate('kinetic_ions/orbits')
kinetic_ions/e1_v1_density/delta_f particle distribution ($\delta f$) out.evaluate('kinetic_ions/delta_f', dataset='kinetic_ions/e1_v1_density/delta_f')
kinetic_ions/e1_v1_density/f particle distribution ($f$) out.evaluate('kinetic_ions/f', dataset='kinetic_ions/e1_v1_density/f')
total_energy scalar time series (total energy) out.evaluate('scalars', variables='total_energy')
Hints
-----
- t=-1 (index), t=slice(...) or t=0.5 (time value) selects snapshots; the t dimension is kept.
- Other keyword arguments select named coordinates, e.g. component=0 or marker=[0, 1, 2].
- Fields: pass eta1=, eta2=, eta3= (scalars or 1D arrays) to evaluate on a logical grid;
omitted directions default to 0.5. The result carries physical coordinates X, Y, Z.
- Particles: out.info('species/variable') lists alternative datasets for dataset=.
- Orbits are an xarray.Dataset with one (t, marker) variable per quantity, e.g. orbits.x;
each variable's 'description' attribute says what it is.
- Results are xarray objects: use .sel/.isel, .plot(x='X'), or .values for NumPy.
[6]:
print(out.kinetic_ions)
ProductNamespace('kinetic_ions', products=('e1_v1_density/delta_f', 'e1_v1_density/f', 'orbits'))
Reconstructed setup and initial conditions#
out.info() prints every evaluable key with a short description and the exact out.evaluate(...) call that loads it, followed by a few hints on selecting times and coordinates. The run configuration is not part of info(): in out.metadata, serialized initial conditions live on each variable under model → species → variables → initial_conditions. out.initial_conditions gives reconstructed objects, including saved Python functions and classes; out.model restores them onto
the model variables. Unsupported definitions remain available in out.metadata.
[7]:
print("Model parameters:", out.model.params)
print("Kinetic variables:", out.model.kinetic_ions.variables)
print("Propagator options:")
for name, options in out.metadata["model"]["propagator_options"].items():
print(f" {name}: {options}")
saved = out.metadata["model"]["species"]["kinetic_ions"]["variables"]["var"]["initial_conditions"]
print("Initial-condition entries:", tuple(saved))
initial = out.initial_conditions["kinetic_ions"]["var"]
print("Background distribution:", initial["backgrounds"])
print("Initial distribution:", initial["initial_condition"])
Model parameters: {'base_units': BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None), 'charge_number': 1, 'mass_number': 1.0, 'alpha': 1.0, 'epsilon': -1.0, 'with_B0': False}
Kinetic variables: {'var': PICVariable (Particles6D)}
Propagator options:
push_eta: {'butcher': {'algo': 'rk4'}}
coupling_va: {'solver': 'pcg', 'precond': 'MassMatrixPreconditioner', 'solver_params': {'tol': 1e-08, 'maxiter': 3000, 'info': False, 'recycle': True}}
Initial-condition entries: ('backgrounds', 'perturbations', 'initial_condition')
Background distribution: Maxwellian3D(
n=(1.0, None),
u1=(0.0, None),
u2=(0.0, None),
u3=(0.0, None),
vth1=(1.0, None),
vth2=(1.0, None),
vth3=(1.0, None),
uniform_on_disc=False,
)
Initial distribution: Maxwellian3D(
n=(1.0, ModesCos(
ls=(1,),
ms=None,
ns=None,
amps=(0.001,),
Lx=1.0,
Ly=1.0,
Lz=1.0,
given_in_basis=None,
comp=0,
perb_domain=(None, None, None),
)),
u1=(0.0, None),
u2=(0.0, None),
u3=(0.0, None),
vth1=(1.0, None),
vth2=(1.0, None),
vth3=(1.0, None),
uniform_on_disc=False,
)
The reconstructed kinetic distributions are ordinary Struphy background objects. Here we evaluate the saved equilibrium and initial distribution along \(\eta_1\); their difference is the density perturbation that seeded the run.
[8]:
background = initial["backgrounds"]
initial_distribution = initial["initial_condition"]
eta1 = np.linspace(0.0, 1.0, 256)
zeros = np.zeros_like(eta1)
n_background = np.asarray(background.n(eta1, zeros, zeros))
n_initial = np.asarray(initial_distribution.n(eta1, zeros, zeros))
fig, ax = plt.subplots()
ax.plot(eta1, n_background, label="background density")
ax.plot(eta1, n_initial, label="initial density")
ax.plot(eta1, n_initial - n_background, label="density perturbation")
ax.set(xlabel=r"$\eta_1$", ylabel="density", title="Saved kinetic initial condition")
ax.legend();
[9]:
phase_space = out.kinetic_ions.e1_v1_density.f
print(phase_space)
<xarray.DataArray 'f' (t: 101, eta1: 32, v1: 32)> Size: 827kB
[103424 values with dtype=float64]
Coordinates:
* t (t) float64 808B 0.0 0.05 0.1 0.15 0.2 ... 4.8 4.85 4.9 4.95 5.0
* eta1 (eta1) float64 256B 0.01562 0.04688 0.07812 ... 0.9531 0.9844
* v1 (v1) float64 256B -4.844 -4.531 -4.219 ... 4.219 4.531 4.844
t_seconds (t) float64 808B 0.0 1.668e-10 3.336e-10 ... 1.651e-08 1.668e-08
Attributes:
label: $f$
long_name: $f$
run: dt=0.05, algo=LieTrotter, Nel=(16, 1, 1), p=(2, 1, 1)
run_name: vlasov_ampere_demo
Evaluate saved splines directly: 1-D, 2-D, and 3-D#
evaluate() reads saved FEEC coefficients and evaluates the spline directly, without creating a post-processing field. With no eta arguments, it evaluates the full simulation grid at its cell centres. For cuts, pass a scalar for each direction to hold fixed; omitted directions stay on the cell centres of the simulation grid. Scalars, lists, NumPy arrays, and range objects can be mixed; every non-scalar input becomes a dimension of the tensor-product result. Thus one varying coordinate
makes a 1-D line, two make a 2-D plane, and no coordinates makes a 3-D volume. The result is always an xarray array with logical eta1/eta2/eta3 coordinates and mapped physical X/Y/Z coordinates.
A representation conversion is applied after spline evaluation. Its input is inferred from the saved FEEC space, so representation= specifies only the target: "0", "1", "2", "3", "v", or "norm". Scalars default to "0"; vectors default to "norm".
[10]:
# 1-D: a field line at eta2 = eta3 = 0.5.
eta1_line = np.linspace(0.0, 1.0, 128)
phi_line = out.evaluate(
"em_fields/phi",
eta1=eta1_line,
eta2=0.5,
eta3=0.5,
t=-1,
)
phi_line.plot()
print(phi_line.dims, phi_line.shape)
('t', 'eta1') (1, 128)
For a 2-D plane, vary two coordinates and hold the third fixed. This is useful for a cross-section of a 3-D field even when the simulation was run on a coarser grid: the spline is evaluated at the requested points. For a 3-D volume, vary all three coordinates. Keep volume grids modest, then take a plane or line from the labeled result for plotting or further analysis.
[11]:
# 2-D: a logical eta1--eta2 plane at eta3 = 0.5.
phi_plane = out.evaluate(
"em_fields/phi",
eta1=np.linspace(0.0, 1.0, 64),
eta2=np.linspace(0.0, 1.0, 48),
eta3=0.5,
t=-1,
)
phi_plane.plot(x="eta1", y="eta2")
# 3-D: vary all three coordinates; keep volume grids modest.
phi_volume = out.evaluate(
"em_fields/phi",
eta1=np.linspace(0.0, 1.0, 32),
eta2=np.linspace(0.0, 1.0, 16),
eta3=np.linspace(0.0, 1.0, 8),
t=-1,
)
print(phi_volume.dims, phi_volume.shape)
# xarray plots a 2-D slice of the volume; choose the mid-plane by coordinate index.
phi_volume.isel(eta3=phi_volume.sizes["eta3"] // 2).plot(x="eta1", y="eta2")
# Without eta arguments the spline is evaluated at the cell centres of the simulation grid,
# here 16 x 1 x 1 cells, so this run's grid gives a line along eta1.
phi_grid = out.evaluate("em_fields/phi", t=-1)
print(phi_grid.dims, phi_grid.shape)
('t', 'eta1', 'eta2', 'eta3') (1, 32, 16, 8)
('t', 'eta1', 'eta2', 'eta3') (1, 16, 1, 1)
Inspecting out itself#
Most of what out exposes is generated on demand (__getattr__, cached properties), so vars(out) only shows a handful of private cache slots, not the products or methods. Use dir(out) for the flat list IPython’s own tab-completion relies on, and out.info() for the full product tree. The same applies one level down: out.model is not a generic stub but the concrete model class of the run (VlasovAmpereOneSpecies here), so dir(out.model) and print(out.model) already
show that model’s own parameters directly — no separate OutputVlasovAmpereOneSpecies-style class is needed.
[12]:
print("vars(out):", vars(out)) # only private cache slots
vars(out): {'path_out': PosixPath('/tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_demo'), '_time_units': 'normalized', 'comm': <feectools.ddm.mpi.MockComm object at 0x7f89bd8d2500>, '_time': None, '_grids_log': None, '_grids_phy': None, '_scalars': <xarray.Dataset> Size: 4kB
Dimensions: (t: 101)
Coordinates:
* t (t) float64 808B 0.0 0.05 0.1 0.15 ... 4.85 4.9 4.95 5.0
t_seconds (t) float64 808B 0.0 1.668e-10 ... 1.651e-08 1.668e-08
Data variables:
electric_energy (t) float64 808B 1.662e-06 1.652e-06 ... 1.489e-07
kinetic_energy (t) float64 808B 9.661 0.0003104 ... -0.000331 -0.0003375
total_energy (t) float64 808B 9.661 0.0003121 ... -0.0003308 -0.0003374, '_products': {'fields': <struphy.post_processing.output.ProductMapping object at 0x7f89bb7e1090>, 'distributions': <struphy.post_processing.output.ProductMapping object at 0x7f89ba4c3730>, 'densities': <struphy.post_processing.output.ProductMapping object at 0x7f89ba4c34f0>, 'orbits': <struphy.post_processing.output.ProductMapping object at 0x7f89ba4c3700>}, '_label': 'dt=0.05, algo=LieTrotter, Nel=(16, 1, 1), p=(2, 1, 1)', '_profile': None, '_tree': <xarray.DataTree>
Group: /
│ Attributes:
│ schema_version: 2
│ options: {"step": 1, "celldivide": [1, 1, 1], "physical": true, "...
├── Group: /em_fields
│ Dimensions: (t: 101, component: 3, eta1: 17, eta2: 2, eta3: 2)
│ Coordinates:
│ * t (t) float64 808B 0.0 0.05 0.1 0.15 0.2 ... 4.85 4.9 4.95 5.0
│ * component (component) int64 24B 0 1 2
│ * eta1 (eta1) float64 136B 0.0 0.0625 0.125 ... 0.875 0.9375 1.0
│ * eta2 (eta2) float64 16B 0.0 1.0
│ * eta3 (eta3) float64 16B 0.0 1.0
│ X (eta1, eta2, eta3) float64 544B ...
│ Y (eta1, eta2, eta3) float64 544B ...
│ Z (eta1, eta2, eta3) float64 544B ...
│ Data variables:
│ e_field (t, component, eta1, eta2, eta3) float64 165kB ...
│ e_field_xyz (t, component, eta1, eta2, eta3) float64 165kB ...
│ phi (t, eta1, eta2, eta3) float64 55kB ...
│ phi_xyz (t, eta1, eta2, eta3) float64 55kB ...
└── Group: /kinetic_ions
├── Group: /kinetic_ions/orbits
│ Dimensions: (t: 101, marker: 12)
│ Coordinates:
│ * t (t) float64 808B 0.0 0.05 0.1 0.15 0.2 ... 4.8 4.85 4.9 4.95 5.0
│ * marker (marker) int64 96B 0 1 2 3 4 5 6 7 8 9 10 11
│ Data variables:
│ x (t, marker) float64 10kB ...
│ y (t, marker) float64 10kB ...
│ z (t, marker) float64 10kB ...
│ v1 (t, marker) float64 10kB ...
│ v2 (t, marker) float64 10kB ...
│ v3 (t, marker) float64 10kB ...
│ weight (t, marker) float64 10kB ...
│ Attributes:
│ product: orbits
│ label: marker orbits
└── Group: /kinetic_ions/e1_v1_density
Dimensions: (t: 101, eta1: 32, v1: 32)
Coordinates:
* t (t) float64 808B 0.0 0.05 0.1 0.15 0.2 ... 4.8 4.85 4.9 4.95 5.0
* eta1 (eta1) float64 256B 0.01562 0.04688 0.07812 ... 0.9531 0.9844
* v1 (v1) float64 256B -4.844 -4.531 -4.219 -3.906 ... 4.219 4.531 4.844
Data variables:
f (t, eta1, v1) float64 827kB ...
delta_f (t, eta1, v1) float64 827kB ..., '_species': {'em_fields', 'kinetic_ions'}, '_seconds': 3.3356409519815204e-09, '_spline_derham': <struphy.feec.psydac_derham.Derham object at 0x7f89b8180d90>, '_spline_snapshots': {100: {'em_fields': {'e_field': <struphy.feec.psydac_derham.SplineFunction object at 0x7f89b8068bb0>, 'phi': <struphy.feec.psydac_derham.SplineFunction object at 0x7f89b806a650>}}}, 'metadata': {'name': '', 'description': '', 'model': {'model': 'VlasovAmpereOneSpecies', 'params': {'base_units': {'BaseUnits': {'x': 1.0, 'B': 1.0, 'n': 1.0, 'kBT': None}}, 'charge_number': 1, 'mass_number': 1.0, 'alpha': 1.0, 'epsilon': -1.0, 'with_B0': False}, 'species': {'em_fields': {'class': 'EMFields', 'charge_number': 0, 'mass_number': 0, 'alpha': None, 'epsilon': None, 'kappa': None, 'variables': {'e_field': {'class': 'FEECVariable', 'space': 'Hcurl', 'save_data': True, 'initial_conditions': {'backgrounds': None, 'perturbations': None}}, 'phi': {'class': 'FEECVariable', 'space': 'H1', 'save_data': True, 'initial_conditions': {'backgrounds': None, 'perturbations': None}}}}, 'kinetic_ions': {'class': 'KineticIons', 'charge_number': 1, 'mass_number': 1.0, 'alpha': 1.0, 'epsilon': -1.0, 'kappa': None, 'variables': {'var': {'class': 'PICVariable', 'space': 'Particles6D', 'save_data': True, 'n_as_volume_form': False, 'initial_conditions': {'backgrounds': {'type': 'Maxwellian3D', 'params': {'n': [1.0, None], 'u1': [0.0, None], 'u2': [0.0, None], 'u3': [0.0, None], 'vth1': [1.0, None], 'vth2': [1.0, None], 'vth3': [1.0, None], 'uniform_on_disc': False}}, 'perturbations': None, 'initial_condition': {'type': 'Maxwellian3D', 'params': {'n': [1.0, {'type': 'ModesCos', 'params': {'ls': [1], 'ms': None, 'ns': None, 'amps': [0.001], 'Lx': 1.0, 'Ly': 1.0, 'Lz': 1.0, 'given_in_basis': None, 'comp': 0, 'perb_domain': [None, None, None]}}], 'u1': [0.0, None], 'u2': [0.0, None], 'u3': [0.0, None], 'vth1': [1.0, None], 'vth2': [1.0, None], 'vth3': [1.0, None], 'uniform_on_disc': False}}}}}, 'loading_params': {'Np': 5000, 'ppc': 32, 'ppb': None, 'loading': 'pseudo_random', 'seed': 1234, 'moments': [0.0, 0.0, 0.0, 1.0, 1.0, 1.0], 'B0': 2.0, 'spatial': 'uniform', 'specific_markers': None, 'set_zero_velocity': [False, False, False], 'n_quad': 1, 'dir_exrernal': None, 'dir_particles': None, 'dir_particles_abs': None, 'restart_key': None}, 'weights_params': {'control_variate': True, 'reject_weights': False, 'threshold': 0.0}, 'boundary_params': {'bc': ['periodic', 'periodic', 'periodic'], 'bc_refill': None, 'bc_sph': ['periodic', 'periodic', 'periodic'], 'mean_velocity_index': None}, 'sorting_params': {'do_sort': True, 'sorting_frequency': 0, 'boxes_per_dim': [4, 1, 1], 'box_bufsize': 2.0, 'dims_mask': [True, True, True]}, 'saving_params': {'n_markers': 12, 'binning_plots': [{'slice': 'e1_v1', 'n_bins': [32, 32], 'ranges': [[0.0, 1.0], [-5.0, 5.0]], 'divide_by_jac': True, 'output_quantity': 'density'}], 'kernel_density_plots': []}, 'bufsize': 1.0}}, 'propagator_options': {'push_eta': {'butcher': {'algo': 'rk4'}}, 'coupling_va': {'solver': 'pcg', 'precond': 'MassMatrixPreconditioner', 'solver_params': {'tol': 1e-08, 'maxiter': 3000, 'info': False, 'recycle': True}}}, 'initial_conditions_schema_version': 1}, 'params_path': None, 'env': {'out_folders': '/tmp/struphy_postprocessing_jqvufi52', 'sim_folder': 'vlasov_ampere_demo', 'sim_label': None, 'restart': False, 'max_runtime': 300, 'save_step': 1, 'save_restart': False, 'sort_step': 0, 'num_clones': 1}, 'time_opts': {'dt': 0.05, 'Tend': 5.0, 'split_algo': 'LieTrotter'}, 'domain': {'type': 'Cuboid', 'params': {'l1': 0.0, 'r1': 6.283185307179586, 'l2': 0.0, 'r2': 1.0, 'l3': 0.0, 'r3': 1.0}}, 'equil': {'type': 'HomogenSlab', 'params': {'B0x': 0.0, 'B0y': 0.0, 'B0z': 1.0, 'beta': 0.1, 'n0': 1.0}}, 'grid': {'num_elements': [16, 1, 1], 'mpi_dims_mask': [True, True, True]}, 'derham_opts': {'degree': [2, 1, 1], 'bcs': [None, None, None], 'nquads': None, 'nquads_proj': None, 'polar_splines': False, 'local_projectors': False}, 'profiling_opts': {'file_path': None, 'label': None, 'use_likwid': None, 'perf_events': None, 'use_line_profiler': None, 'use_memray': None, 'memory_profile_path': None, 'memray_native_traces': None, 'memray_trace_python_allocators': None, 'memray_follow_fork': None, 'deactivate_profiling': None, 'use_nvtx': None, 'use_gpu_timing': None, 'gpu_timing_backend': None, 'deactivate_file_output': None, 'recursive_profile': None, 'aggregation_mode': None, 'profile_mpi_calls': None, 'track_threads': None, 'track_async': None, 'capture_region_source': None, 'metadata_detail': None, 'buffer_limit': None, 'output_mode': None, 'hdf5_compression': None, 'hdf5_compression_level': None, 'hdf5_chunk_size': None, 'memray': None, 'gpu': None, 'hdf5': None}, 'mpi_ranks': 1, 'use_mpi_comm_world': False, 'started_at_epoch_s': 1790926356.3625028, 'one_time_step': False, 'profiling_activated': True}, 'mpi_ranks': 1, 'grid': TensorProductGrid(
num_elements=(16, 1, 1),
mpi_dims_mask=(True, True, True),
), 'derham_opts': DerhamOptions(
degree=(2, 1, 1),
bcs=(None, None, None),
nquads=None,
nquads_proj=None,
polar_splines=False,
local_projectors=False,
), 'domain': Cuboid(
l1=0.0,
r1=6.283185307179586,
l2=0.0,
r2=1.0,
l3=0.0,
r3=1.0,
), 'initial_conditions': {'em_fields': {'e_field': {'backgrounds': None, 'perturbations': None}, 'phi': {'backgrounds': None, 'perturbations': None}}, 'kinetic_ions': {'var': {'backgrounds': Maxwellian3D(
n=(1.0, None),
u1=(0.0, None),
u2=(0.0, None),
u3=(0.0, None),
vth1=(1.0, None),
vth2=(1.0, None),
vth3=(1.0, None),
uniform_on_disc=False,
), 'perturbations': None, 'initial_condition': Maxwellian3D(
n=(1.0, ModesCos(
ls=(1,),
ms=None,
ns=None,
amps=(0.001,),
Lx=1.0,
Ly=1.0,
Lz=1.0,
given_in_basis=None,
comp=0,
perb_domain=(None, None, None),
)),
u1=(0.0, None),
u2=(0.0, None),
u3=(0.0, None),
vth1=(1.0, None),
vth2=(1.0, None),
vth3=(1.0, None),
uniform_on_disc=False,
)}}}, 'model': VlasovAmpereOneSpecies(
base_units=BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None),
charge_number=1,
mass_number=1.0,
alpha=1.0,
epsilon=-1.0,
with_B0=False,
), 'time_opts': Time(
dt=0.05,
Tend=5.0,
split_algo='LieTrotter',
)}
[13]:
print("dir(out):", [name for name in dir(out) if not name.startswith("_")])
dir(out): ['catalog', 'clear_cache', 'close', 'comm', 'compare', 'densities', 'density_catalog', 'derham_opts', 'distribution_catalog', 'distributions', 'domain', 'em_fields', 'equil', 'evaluate', 'field_catalog', 'fields', 'grid', 'grids_log', 'grids_phy', 'info', 'initial_conditions', 'is_processed', 'iter_spline_coefficients', 'keys', 'kinetic_ions', 'label', 'metadata', 'model', 'mpi_ranks', 'orbit_catalog', 'orbits', 'path_out', 'path_pproc', 'pproc', 'profile', 'provenance', 'report', 'save_scalars', 'scalars', 'seconds_per_time', 'species_catalog', 'spline_fields', 'time', 'time_opts', 'time_scale', 'time_unit', 'time_units', 'to_si', 'tree', 'units', 'with_physical_coords', 'with_time_units', 'xarray']
[14]:
print(out.model) # the concrete model class, with its own parameters
VlasovAmpereOneSpecies(
base_units=BaseUnits(x=1.0, B=1.0, n=1.0, kBT=None),
charge_number=1,
mass_number=1.0,
alpha=1.0,
epsilon=-1.0,
with_B0=False,
)
[15]:
print("dir(out.model):", [name for name in dir(out.model) if not name.startswith("_")][:15])
dir(out.model): ['EMFields', 'KineticIons', 'Propagators', 'base_units', 'bulk_species', 'cannot_be_used_for', 'cannot_be_used_for_html', 'cannot_be_used_for_latex', 'cannot_be_used_for_markdown', 'clone_config', 'create_doc', 'diagnostic_species', 'discretization', 'discretization_html', 'discretization_latex']
Products are xarray arrays#
A product is an xarray.DataArray, so xarray’s own plotting already draws it, with the labels and units Struphy stored:
[16]:
phase_space.isel(t=-1).plot(x="eta1", y="v1")
[16]:
<matplotlib.collections.QuadMesh at 0x7f89b3f37790>
Use xarray’s plotting methods for ordinary one- and two-dimensional output. Select a saved time or spatial slice with .isel() or .sel() first, then call .plot() or .plot.line().
Scalar overview and time series#
Field and particle products are xarray DataArray objects; out.evaluate("scalars") returns an xarray Dataset. out.evaluate("kinetic_ions/f", dataset="e1_v1_density/f") looks a product up by name, which suits scripts and loops; attribute access is convenient interactively.
Call .plot() for a scalar time series. Use a Matplotlib axes when combining several series or setting plot options.
[17]:
out.scalars.electric_energy.plot.line(x="t")
[17]:
[<matplotlib.lines.Line2D at 0x7f89b3e199f0>]
[18]:
t_fit = 2.0 # Struphy time units, like every time coordinate of this run
energy = out.scalars.electric_energy.sel(t=slice(0.0, t_fit))
fig, ax = plt.subplots()
energy.plot.line(ax=ax, label="electric energy")
ax.set_yscale("log")
ax.legend()
[18]:
<matplotlib.legend.Legend at 0x7f89b3e3b0d0>
Two-dimensional data#
Choose the displayed dimensions with x and y, and pick one value for every other dimension by naming it: t="last" (or "first"), t=-1 for a position, and t=0.35 for the nearest coordinate value. Arrays can also be sliced beforehand with xarray’s .isel() and .sel(). coords="physical" draws on the mapped coordinates instead of the logical ones.
[19]:
phase_space.isel(t=-1).plot(x="eta1", y="v1")
[19]:
<matplotlib.collections.QuadMesh at 0x7f89b3d8ff70>
For a compact view of the evolution, select saved times and use xarray faceting.
[20]:
phase_space.isel(t=np.linspace(0, phase_space.sizes["t"] - 1, 5, dtype=int)).plot(
x="eta1", y="v1", col="t", col_wrap=5
)
[20]:
<xarray.plot.facetgrid.FacetGrid at 0x7f89b3fbf250>
Selecting saved snapshots#
Use .isel() for index-based selection and .sel() for coordinate-based selection. This keeps selection explicit and works with every xarray operation.
[21]:
final_phase_space = phase_space.isel(t=-1)
final_phase_space.plot(x="eta1", y="v1")
[21]:
<matplotlib.collections.QuadMesh at 0x7f89b3910220>
Saved marker orbits sit under their species as an xarray.Dataset with one (t, marker) variable per quantity (x, y, z, velocities, weight); each variable’s description attribute says what it is. Select one marker and plot its positions over time.
[22]:
orbit = out.kinetic_ions.orbits.isel(marker=0)
orbit[["x", "y", "z"]].to_dataarray("quantity").plot.line(x="t", hue="quantity")
[22]:
[<matplotlib.lines.Line2D at 0x7f89b3971f30>,
<matplotlib.lines.Line2D at 0x7f89b3971000>,
<matplotlib.lines.Line2D at 0x7f89b3970e80>]
A small collection of explicit snapshots is often more useful in a reproducible notebook than an interactive widget or animation.
[23]:
phase_space.isel(t=[0, -1]).plot(x="eta1", y="v1", col="t")
[23]:
<xarray.plot.facetgrid.FacetGrid at 0x7f89b3f82cb0>
[24]:
fig, ax = plt.subplots()
phase_space.isel(t=-1).plot(ax=ax, x="eta1", y="v1")
fig.savefig(os.path.join(demo_root, "phase_space_final.png"), bbox_inches="tight")
The reconstructed equilibrium is available directly on the output handle for inspection and for model-specific analysis.
[25]:
print(out.equil)
HomogenSlab(
B0x=0.0,
B0y=0.0,
B0z=1.0,
beta=0.1,
n0=1.0,
)
Derived quantities#
xarray arithmetic computes derived quantities without special APIs. For example, subtract the first sample to obtain a drift and divide by it to obtain a relative error.
[26]:
total_energy = out.scalars.total_energy
energy_drift = total_energy - total_energy.isel(t=0)
energy_error = abs(energy_drift) / abs(total_energy.isel(t=0))
print(f"largest drift of the total energy: {abs(energy_drift).max().item():.3e}")
energy_error.plot.line(x="t")
largest drift of the total energy: 9.661e+00
[26]:
[<matplotlib.lines.Line2D at 0x7f89bb7cce80>]
For custom diagnostics, first select the labeled subset needed for the calculation. Here a space-time field line is retained as an xarray object, ready for NumPy, SciPy, or another analysis package.
[27]:
space_time_line = out.evaluate(
"em_fields/phi", eta1=np.linspace(0.0, 1.0, 64), eta2=0.5, eta3=0.5
)
print(space_time_line.dims, space_time_line.shape)
('t', 'eta1') (101, 64)
Reducing distribution functions#
A binned distribution usually has more dimensions than the question needs. Xarray reductions retain the remaining named dimensions, so averaging over e1, e2 and e3 turns the \((\eta_1, v_1)\) product into \(f(v_1, t)\). The mean is uniform in logical coordinates, which is a volume average on a Cartesian domain.
[28]:
# average over whichever logical space directions the binning kept
spatial = [dim for dim in ("eta1", "eta2", "eta3") if dim in phase_space.dims]
f_of_v = phase_space.mean(spatial)
print(f_of_v.dims)
f_of_v.plot(x="t", y="v1")
('t', 'v1')
[28]:
<matplotlib.collections.QuadMesh at 0x7f89b3722fb0>
Velocity moments are weighted xarray reductions. The bin widths and velocity coordinate remain labeled, making the density, mean velocity, and variance explicit.
[29]:
dv1 = phase_space.v1.differentiate("v1")
density = (phase_space * dv1).sum("v1")
mean_v1 = (phase_space * phase_space.v1 * dv1).sum("v1") / density
variance_v1 = (phase_space * (phase_space.v1 - mean_v1) ** 2 * dv1).sum("v1") / density
mean_density = density.mean(spatial)
mean_density.plot.line(x="t")
variance_v1.mean(spatial).plot.line(x="t")
[29]:
[<matplotlib.lines.Line2D at 0x7f89b37fbbe0>]
Physical units#
Products are in the normalization of the model. out.units holds the units of that normalization and out.to_si(product) converts a product: the time coordinate to seconds, the mapped coordinates X, Y, Z to meters and the velocities v1, v2, v3 to m/s. The values are converted only when unit= names the unit the variable was normalized with ("x", "B", "n", "v", "t", "p", "rho", "j" or "kBT"), or a number for a composite unit
together with its label, because a product does not record which unit its variable uses. The original product is not modified.
[30]:
print(f"1 length unit = {out.units.x} m, 1 velocity unit = {out.units.v:.4g} m/s, 1 time unit = {out.units.t:.4g} s")
phase_space_si = out.to_si(phase_space)
print(phase_space_si.v1.attrs["units"], phase_space_si.t.attrs["units"])
phase_space_si.isel(t=-1).plot(x="eta1", y="v1")
1 length unit = 1.0 m, 1 velocity unit = 2.998e+08 m/s, 1 time unit = 3.336e-09 s
m/s s
[30]:
<matplotlib.collections.QuadMesh at 0x7f89b363c160>
Save standard output#
Products are xarray objects, so figures are saved with Matplotlib’s plt.savefig(path) and arrays with xarray’s own writers, for example phi.to_netcdf(path). For a compact record of the run, out.report() writes the full scalar history as scalars.csv and a report listing the run metadata, every available product, and the dimensions and units of any products passed with products=. By default the files go to post_processing/report/; format is "markdown" or "html",
and the return value is the path of the report file.
[31]:
report = out.report(products=["em_fields/phi"], format="html")
print("Wrote:")
for path in (report, os.path.join(os.path.dirname(report), "scalars.csv")):
print(" ", os.path.relpath(path, out.path_out))
Wrote:
post_processing/report/report.html
post_processing/report/scalars.csv
Comparing runs#
Time series accept arrays of other simulations, so comparing runs needs nothing special. Series are labelled by the run they come from, and the runs may have different time grids.
[32]:
sim_coarse = Simulation(
model=build_model(),
env=EnvironmentOptions(out_folders=demo_root, sim_folder="vlasov_ampere_coarse", save_restart=False),
time_opts=Time(dt=0.1, Tend=5.0),
domain=domains.Cuboid(r1=2 * 3.141592653589793),
equil=equils.HomogenSlab(),
grid=grids.TensorProductGrid(num_elements=(16, 1, 1)),
derham_opts=DerhamOptions(degree=(2, 1, 1)),
)
out_coarse = sim_coarse.run(profiling_activated=True)
fig, ax = plt.subplots()
out.scalars.electric_energy.plot.line(ax=ax, label="dt = 0.05")
out_coarse.scalars.electric_energy.plot.line(ax=ax, label="dt = 0.1")
ax.legend()
/opt/hostedtoolcache/Python/3.10.21/x64/lib/python3.10/site-packages/struphy/models/species.py:215: UserWarning: Override equation parameter self.alpha =1.0
warnings.warn(f"Override equation parameter {self.alpha =}")
/opt/hostedtoolcache/Python/3.10.21/x64/lib/python3.10/site-packages/struphy/models/species.py:222: UserWarning: Override equation parameter self.epsilon =-1.0
warnings.warn(f"Override equation parameter {self.epsilon =}")
Stabilizing Poisson solve with self.options.sigma_1 =1e-14
Time stepping: 100%|██████████| 50/50 [00:00<00:00, 84.67step/s]
╭────────────────────────────────────────────────────────────────────────────────────────────╮
│ region % session total [s] │
├────────────────────────────────────────────────────────────────────────────────────────────┤
│ scope_profiler.session 100.00% 0.745539 │
│ └─ (own) 1.08% 0.008020 │
│ └─ setup: total 20.56% 0.153303 │
│ │ └─ (own) 0.09% 0.000695 │
│ │ └─ setup: allocate 9.98% 0.074435 │
│ │ │ └─ (own) 0.01% 0.000071 │
│ │ │ └─ setup: feec 5.83% 0.043456 │
│ │ │ │ └─ (own) 0.01% 0.000095 │
│ │ │ │ └─ setup: derham 5.81% 0.043315 │
│ │ │ │ └─ setup: mass ops 0.00% 0.000008 │
│ │ │ │ └─ setup: basis ops 0.00% 0.000008 │
│ │ │ │ └─ setup: projected equil 0.00% 0.000029 │
│ │ │ └─ setup: variables 0.39% 0.002916 │
│ │ │ │ └─ (own) 0.01% 0.000077 │
│ │ │ │ └─ setup var: em_fields.e_field 0.02% 0.000130 │
│ │ │ │ └─ setup var: em_fields.phi 0.01% 0.000043 │
│ │ │ │ └─ setup var: kinetic_ions.var 0.36% 0.002666 │
│ │ │ │ │ └─ (own) 0.33% 0.002471 │
│ │ │ │ │ └─ do_sort 0.03% 0.000195 │
│ │ │ │ │ │ └─ (own) 0.01% 0.000042 │
│ │ │ │ │ │ └─ put_particles_in_boxes 0.02% 0.000152 │
│ │ │ └─ setup: propagators 2.72% 0.020313 │
│ │ │ │ └─ (own) 0.01% 0.000039 │
│ │ │ │ └─ setup prop: PushEta 0.01% 0.000038 │
│ │ │ │ └─ setup prop: VlasovAmpereCoupling 2.71% 0.020236 │
│ │ │ └─ setup: helpers 1.03% 0.007680 │
│ │ │ │ └─ (own) 0.45% 0.003377 │
│ │ │ │ └─ accum: charge_density_0form 0.04% 0.000268 │
│ │ │ │ │ └─ (own) 0.00% 0.000029 │
│ │ │ │ │ └─ kernel: charge_density_0form 0.02% 0.000185 │
│ │ │ │ │ └─ accum comm: charge_density_0form 0.01% 0.000054 │
│ │ │ │ └─ solve: PoissonSolve 0.53% 0.003983 │
│ │ │ │ └─ update_feec_variables 0.01% 0.000052 │
│ │ └─ setup: run metadata 0.15% 0.001119 │
│ │ └─ setup: data storage 0.71% 0.005285 │
│ │ └─ setup: geometry vtk 7.84% 0.058417 │
│ │ └─ setup: plasma params 0.39% 0.002931 │
│ │ └─ setup: initial diagnostics 0.13% 0.000973 │
│ │ └─ setup: hdf5 datasets 1.27% 0.009450 │
│ └─ model.integrate (50x) 37.31% 0.278135 │
│ │ └─ (own) 0.08% 0.000594 │
│ │ └─ prop: PushEta (50x) 4.59% 0.034249 │
│ │ │ └─ (own) 2.29% 0.017054 │
│ │ │ └─ pusher: push_eta_stage (50x) 2.31% 0.017194 │
│ │ │ │ └─ (own) 0.97% 0.007230 │
│ │ │ │ └─ kernel: push_eta_stage (200x) 1.34% 0.009964 │
│ │ └─ prop: VlasovAmpereCoupling (50x) 32.63% 0.243293 │
│ │ │ └─ (own) 2.84% 0.021164 │
│ │ │ └─ accum: vlasov_maxwell (50x) 9.78% 0.072939 │
│ │ │ │ └─ (own) 0.46% 0.003448 │
│ │ │ │ └─ kernel: vlasov_maxwell (50x) 4.54% 0.033853 │
│ │ │ │ └─ accum comm: vlasov_maxwell (100x) 4.78% 0.035639 │
│ │ │ └─ solve: SchurSolver (50x) 17.52% 0.130592 │
│ │ │ └─ pusher: push_v_with_efield (50x) 1.73% 0.012879 │
│ │ │ │ └─ (own) 0.51% 0.003834 │
│ │ │ │ └─ kernel: push_v_with_efield (50x) 1.21% 0.009045 │
│ │ │ └─ update_feec_variables (50x) 0.77% 0.005718 │
│ └─ diagnostics (50x) 4.97% 0.037051 │
│ └─ save data (51x) 36.09% 0.269030 │
╰────────────────────────────────────────────────────────────────────────────────────────────╯
╭─ Info ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────╮
│ Summary: ../../../../../../../../tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_coarse/profiling_data.h5 (1 rank) │
│ │
│ Explore: │
│ Inspect: scope-profiler inspect ../../../../../../../../tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_coarse/profiling_data.h5 │
│ TUI: scope-profiler tui ../../../../../../../../tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_coarse/profiling_data.h5 │
│ │
│ Visualize and export: │
│ Plot: scope-profiler plot default ../../../../../../../../tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_coarse/profiling_data.h5 -o plots --show │
│ Report: scope-profiler report ../../../../../../../../tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_coarse/profiling_data.h5 -o report.html │
│ Export: scope-profiler export plot-data ../../../../../../../../tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_coarse/profiling_data.h5 -o data │
│ Lines: scope-profiler line-profile ../../../../../../../../tmp/struphy_postprocessing_jqvufi52/vlasov_ampere_coarse/profiling_data.h5 │
│ │
│ Compare runs: │
│ Diff: scope-profiler diff BASE.h5 CANDIDATE.h5 │
│ Check: scope-profiler check BASE.h5 CANDIDATE.h5 │
│ │
│ Durations are in seconds. │
│ Regions may nest, so the summed total can exceed the wall-clock time. │
│ % session uses wall-clock coverage; overlapping recursive calls count once. │
╰───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────╯
[32]:
<matplotlib.legend.Legend at 0x7f89b340c9d0>
Profiling#
A run started with sim.run(profiling_activated=True) records how long each region of the code takes, and out.profile reads that record back. Regions are the setup steps, every propagator (prop: ...), pusher, accumulation, compiled kernel (kernel: ...) and linear solve. Regions nest, so a region’s time includes the regions it calls and the times of different regions must not be added up. summary() returns a dataset along the dimension region, sorted by total time;
table() prints it. Filter with prefix and limit with top.
[33]:
print(out.profile.table(top=8))
kernels = out.profile.summary(prefix="kernel:")
print(kernels.total_time.to_series())
Profile: dt=0.05, algo=LieTrotter, Nel=(16, 1, 1), p=(2, 1, 1) (1 rank(s), 1.319 s)
Region Calls Total [s] Mean [ms] Share
-------------------------- -------- ---------- ---------- -------
scope_profiler.session 1 1.319 1319.101 100.0%
model.integrate 100 0.548 5.485 41.6%
save data 101 0.529 5.239 40.1%
prop: VlasovAmpereCoupling 100 0.482 4.819 36.5%
solve: SchurSolver 100 0.258 2.578 19.5%
setup: total 1 0.154 154.184 11.7%
accum: vlasov_maxwell 100 0.144 1.444 10.9%
setup: allocate 1 0.077 76.512 5.8%
region
kernel: vlasov_maxwell 0.067735
kernel: push_eta_stage 0.019063
kernel: push_v_with_efield 0.018398
kernel: charge_density_0form 0.000142
Name: total_time, dtype: float64
compare() puts the same statistic of several runs side by side, with runs whose region is missing as NaN. Here the two runs of the previous section differ only in the time step, so the number of calls per propagator halves for dt = 0.1.
[34]:
calls = out.profile.compare(out_coarse, metric="calls", prefix="prop:")
print(calls)
<xarray.DataArray 'calls' (run: 2, region: 2)> Size: 32B
array([[100., 100.],
[ 50., 50.]])
Coordinates:
* region (region) <U26 208B 'prop: VlasovAmpereCoupling' 'prop: PushEta'
* run (run) <U53 424B 'dt=0.05, algo=LieTrotter, Nel=(16, 1, 1), p=(2,...
Attributes:
label: calls per rank
Other models#
The interface is the same for every model; only the products differ. Two more short runs show the two product types the Vlasov–Ampère demo does not have: SPH densities, and vector fields on a mapped domain.
SPH densities#
A standing sound wave discretized with SPH markers. KernelDensityPlot reconstructs the density on a grid, which appears under out.densities, while BinningPlot produces the binned quantities under out.distributions.
[35]:
sph_model = ViscousEulerSPH(with_B0=False, with_viscosity=False)
sph_model.propagators.push_eta.options = sph_model.propagators.push_eta.Options(
butcher=ButcherTableau(algo="forward_euler"),
)
sph_model.propagators.push_sph_p.options = sph_model.propagators.push_sph_p.Options(kernel_type="gaussian_1d")
sph_model.euler_fluid.set_markers(
loading_params=LoadingParameters(ppb=8, loading="tesselation"),
weights_params=WeightsParameters(),
boundary_params=BoundaryParameters(),
sorting_params=SortingParameters(boxes_per_dim=(12, 1, 1), dims_mask=(True, False, False)),
saving_params=SavingParameters(
binning_plots=(BinningPlot(slice="e1", n_bins=(32,), ranges=(0.0, 1.0)),),
kernel_density_plots=(KernelDensityPlot(pts_e1=41, pts_e2=1),),
),
)
sph_model.euler_fluid.var.add_background(equils.ConstantVelocity())
sph_model.euler_fluid.var.add_perturbation(del_n=perturbations.ModesSin(ls=(1,), amps=(1.0e-2,)))
sph = Simulation(
model=sph_model,
env=EnvironmentOptions(out_folders=demo_root, sim_folder="sph_soundwave", save_restart=False),
time_opts=Time(dt=0.03125, Tend=2.5, split_algo="Strang"),
domain=domains.Cuboid(r1=2.5),
grid=None,
derham_opts=None,
)
out_sph = sph.run()
print("densities:", tuple(out_sph.density_catalog))
print("binned:", tuple(out_sph.distribution_catalog))
Time stepping: 100%|██████████| 80/80 [00:00<00:00, 147.46step/s]
No post-processed data in /tmp/struphy_postprocessing_jqvufi52/sph_soundwave, processing with default options (call out.pproc(...) to choose them)
Post-processing path /tmp/struphy_postprocessing_jqvufi52/sph_soundwave
No feec fields found in hdf5 file, skipping post-processing of fields.
Evaluation of 3 marker orbits for euler_fluid
100%|██████████| 81/81 [00:00<00:00, 912.10it/s]
Evaluation of distribution functions for euler_fluid
0 starting post-processing of distribution functions for /tmp/struphy_postprocessing_jqvufi52/sph_soundwave/post_processing/kinetic_data/euler_fluid ...
100%|██████████| 1/1 [00:00<00:00, 2492.16it/s]
0%| | 0/1 [00:00<?, ?it/s]rank = 0 ----------------------------
self._pproc_rank =0 with xp.sum(data) =np.float64(1036.8000000000002) and xp.sum(data_df) =np.float64(1036.8000000000002)
self._pproc_rank =0 with xp.sum(data) =np.float64(1036.8000000000002) and xp.sum(data_df) =np.float64(1036.8000000000002)
self._pproc_rank =0 done.
100%|██████████| 1/1 [00:00<00:00, 41.90it/s]
Evaluation of sph density for euler_fluid
100%|██████████| 1/1 [00:00<00:00, 23.83it/s]
densities: ('euler_fluid/view_0/n',)
binned: ('euler_fluid/e1_density/f', 'euler_fluid/e1_density/delta_f')
For a one-dimensional run, the clearest picture is a space-time map: the sweep dimension t may be used as a display axis.
[36]:
density = out_sph.euler_fluid.view_0.n.isel(eta2=0, eta3=0)
density.plot(x="t", y="eta1")
[36]:
<matplotlib.collections.QuadMesh at 0x7f89b33df640>
Products are plain xarray.DataArray objects, so anything xarray can do works directly, for example profiles at selected times:
[37]:
density.isel(t=[0, len(density.t) // 4, len(density.t) // 2]).plot.line(x="eta1")
[37]:
[<matplotlib.lines.Line2D at 0x7f89b3991570>,
<matplotlib.lines.Line2D at 0x7f89b32ccac0>,
<matplotlib.lines.Line2D at 0x7f89b32ccc10>]
Vector fields on a mapped domain#
A coaxial waveguide mode of the Maxwell model, on an annulus. With physical=True the post-processing also computes the Cartesian field components (*_xyz), and coords="physical" draws them on the mapped grid, with the plane chosen by plane.
[38]:
a1, a2 = 2.326744, 3.686839
maxwell_model = Maxwell()
maxwell_model.propagators.maxwell.options = maxwell_model.propagators.maxwell.Options(algo="implicit")
maxwell_model.em_fields.e_field.add_perturbation(perturbations.CoaxialWaveguideElectric_r(m=3, a1=a1, a2=a2))
maxwell_model.em_fields.e_field.add_perturbation(perturbations.CoaxialWaveguideElectric_theta(m=3, a1=a1, a2=a2))
maxwell_model.em_fields.b_field.add_perturbation(perturbations.CoaxialWaveguideMagnetic(m=3, a1=a1, a2=a2))
coaxial = Simulation(
model=maxwell_model,
env=EnvironmentOptions(out_folders=demo_root, sim_folder="coaxial", save_restart=False),
time_opts=Time(dt=0.05, Tend=2.0),
domain=domains.HollowCylinder(a1=a1, a2=a2, Lz=2.0),
equil=equils.HomogenSlab(),
grid=grids.TensorProductGrid(num_elements=(24, 48, 1)),
derham_opts=DerhamOptions(degree=(2, 2, 1), bcs=(("dirichlet", "dirichlet"), None, None)),
)
out_coaxial = coaxial.run()
out_coaxial.pproc(physical=True)
print("fields:", tuple(out_coaxial.field_catalog))
print("dimensions:", out_coaxial.em_fields.b_field_xyz.dims)
Time stepping: 100%|██████████| 40/40 [00:00<00:00, 41.34step/s]
Post-processing path /tmp/struphy_postprocessing_jqvufi52/coaxial
Reading hdf5 data of following species:
em_fields:
b_field: <HDF5 group "/feec/em_fields/b_field" (3 members)>
e_field: <HDF5 group "/feec/em_fields/e_field" (3 members)>
Creation of Struphy Fields done.
Evaluating fields ...
100%|██████████| 41/41 [00:00<00:00, 107.91it/s]
No kinetic data found in hdf5 file, skipping post-processing of kinetic data.
fields: ('em_fields/b_field', 'em_fields/b_field_xyz', 'em_fields/e_field', 'em_fields/e_field_xyz')
dimensions: ('t', 'component', 'eta1', 'eta2', 'eta3')
[39]:
out_coaxial.em_fields.b_field_xyz.isel(t=-1, component=2, eta3=0).plot(x="eta1", y="eta2")
[39]:
<matplotlib.collections.QuadMesh at 0x7f89b108eef0>
[40]:
out_coaxial.em_fields.b_field_xyz.isel(component=2, eta3=0).plot(x="eta1", y="eta2", col="t", col_wrap=4)
[40]:
<xarray.plot.facetgrid.FacetGrid at 0x7f89bb7e1fc0>
Representation conversion on a torus#
A toroidal map makes the distinction between FEEC representations visible. The saved electric field is an H(curl) 1-form, so its source representation is inferred as 1. We evaluate one poloidal line directly from saved spline coefficients, then request its native 1-form, normalized-vector, and Cartesian-vector representations. These differ away from a Cartesian map because the metric factors vary around the torus.
[41]:
torus_model = Maxwell()
torus_model.em_fields.e_field.save_data = True
torus_model.em_fields.e_field.add_perturbation(
perturbations.ModesCos(ms=(1,), amps=(0.1,), given_in_basis="1", comp=1)
)
torus = Simulation(
model=torus_model,
env=EnvironmentOptions(out_folders=demo_root, sim_folder="representation_torus", save_restart=False),
time_opts=Time(dt=0.05, Tend=0.05),
domain=domains.HollowTorus(a1=0.2, a2=0.4, R0=1.0, tor_period=1),
equil=equils.HomogenSlab(),
grid=grids.TensorProductGrid(num_elements=(6, 12, 2)),
derham_opts=DerhamOptions(degree=(2, 2, 2), bcs=(("dirichlet", "dirichlet"), None, None)),
)
out_torus = torus.run()
Time stepping: 100%|██████████| 1/1 [00:00<00:00, 16.03step/s]
[42]:
eta2_line = np.linspace(0.0, 1.0, 256)
component=0
common = dict(eta1=0.7, eta2=eta2_line, eta3=0.0, t=-1, component=component)
e_1 = out_torus.evaluate("em_fields/e_field", representation="1", **common)
e_norm = out_torus.evaluate("em_fields/e_field", representation="norm", **common)
e_v = out_torus.evaluate("em_fields/e_field", representation="v", **common)
fig, ax = plt.subplots(ncols=3, figsize=(15, 5))
i = 0
for field, label in ((e_1, "1-form"), (e_norm, "normalized vector"), (e_v, "Vector field")):
print(label, field.isel(t=0))
ax[i].plot(eta2_line, field.isel(t=0), label=label)
ax[i].set(xlabel=r"$\eta_2$", ylabel=label, title="H(curl) field representations on a torus")
i += 1
# ax.legend()
1-form <xarray.DataArray 'e_field' (eta2: 256)> Size: 2kB
array([ 2.51162210e-17, -3.53796800e-06, -7.04305545e-06, -1.05152623e-05,
-1.39545887e-05, -1.73610345e-05, -2.07345997e-05, -2.40752844e-05,
-2.73830886e-05, -3.06580122e-05, -3.39000552e-05, -3.71092177e-05,
-4.02854997e-05, -4.34289011e-05, -4.65394219e-05, -4.96170622e-05,
-5.26618219e-05, -5.56737011e-05, -5.86526998e-05, -6.15988179e-05,
-6.45120554e-05, -6.73924124e-05, -7.02251026e-05, -7.29739818e-05,
-7.56374070e-05, -7.82153784e-05, -8.07078959e-05, -8.31149595e-05,
-8.54365691e-05, -8.76727249e-05, -8.98234267e-05, -9.18886747e-05,
-9.38684687e-05, -9.57628089e-05, -9.75716951e-05, -9.92951274e-05,
-1.00933106e-04, -1.02485630e-04, -1.03952701e-04, -1.05334318e-04,
-1.06630480e-04, -1.07841189e-04, -1.08966444e-04, -1.10004504e-04,
-1.10944919e-04, -1.11785949e-04, -1.12527593e-04, -1.13169851e-04,
-1.13712723e-04, -1.14156209e-04, -1.14500310e-04, -1.14745025e-04,
-1.14890354e-04, -1.14936298e-04, -1.14882856e-04, -1.14730028e-04,
-1.14477814e-04, -1.14126214e-04, -1.13675229e-04, -1.13124858e-04,
-1.12475101e-04, -1.11725958e-04, -1.10877430e-04, -1.09929516e-04,
-1.08883540e-04, -1.07768622e-04, -1.06596676e-04, -1.05367702e-04,
-1.04081700e-04, -1.02738670e-04, -1.01338611e-04, -9.98815236e-05,
-9.83674083e-05, -9.67962647e-05, -9.51680928e-05, -9.34828926e-05,
-9.17406642e-05, -8.99414076e-05, -8.80851226e-05, -8.61718094e-05,
...
8.61718094e-05, 8.80851226e-05, 8.99414076e-05, 9.17406642e-05,
9.34828926e-05, 9.51680928e-05, 9.67962647e-05, 9.83674083e-05,
9.98815236e-05, 1.01338611e-04, 1.02738670e-04, 1.04081700e-04,
1.05367702e-04, 1.06596676e-04, 1.07768622e-04, 1.08883540e-04,
1.09929516e-04, 1.10877430e-04, 1.11725958e-04, 1.12475101e-04,
1.13124858e-04, 1.13675229e-04, 1.14126214e-04, 1.14477814e-04,
1.14730028e-04, 1.14882856e-04, 1.14936298e-04, 1.14890354e-04,
1.14745025e-04, 1.14500310e-04, 1.14156209e-04, 1.13712723e-04,
1.13169851e-04, 1.12527593e-04, 1.11785949e-04, 1.10944919e-04,
1.10004504e-04, 1.08966444e-04, 1.07841189e-04, 1.06630480e-04,
1.05334318e-04, 1.03952701e-04, 1.02485630e-04, 1.00933106e-04,
9.92951274e-05, 9.75716951e-05, 9.57628089e-05, 9.38684687e-05,
9.18886747e-05, 8.98234267e-05, 8.76727249e-05, 8.54365691e-05,
8.31149595e-05, 8.07078959e-05, 7.82153784e-05, 7.56374070e-05,
7.29739818e-05, 7.02251026e-05, 6.73924124e-05, 6.45120554e-05,
6.15988179e-05, 5.86526998e-05, 5.56737011e-05, 5.26618219e-05,
4.96170622e-05, 4.65394219e-05, 4.34289011e-05, 4.02854997e-05,
3.71092177e-05, 3.39000552e-05, 3.06580122e-05, 2.73830886e-05,
2.40752844e-05, 2.07345997e-05, 1.73610345e-05, 1.39545887e-05,
1.05152623e-05, 7.04305545e-06, 3.53796800e-06, 2.51009744e-17])
Coordinates:
t float64 8B 0.05
* eta2 (eta2) float64 2kB 0.0 0.003922 0.007843 ... 0.9922 0.9961 1.0
component int64 8B 0
t_seconds float64 8B 1.668e-10
X (eta2) float64 2kB 1.34 1.34 1.34 1.339 ... 1.339 1.34 1.34 1.34
Y (eta2) float64 2kB -0.0 -0.0 -0.0 -0.0 ... -0.0 -0.0 -0.0 -0.0
Z (eta2) float64 2kB 0.0 0.008377 0.01675 ... -0.008377 -8.328e-17
Attributes:
run: dt=0.05, algo=LieTrotter, Nel=(6, 12, 2), p=(2, 2, 2)
run_name: representation_torus
normalized vector <xarray.DataArray 'e_field' (eta2: 256)> Size: 2kB
array([ 1.25581105e-16, -1.76898400e-05, -3.52152773e-05, -5.25763117e-05,
-6.97729435e-05, -8.68051725e-05, -1.03672999e-04, -1.20376422e-04,
-1.36915443e-04, -1.53290061e-04, -1.69500276e-04, -1.85546089e-04,
-2.01427498e-04, -2.17144505e-04, -2.32697109e-04, -2.48085311e-04,
-2.63309110e-04, -2.78368506e-04, -2.93263499e-04, -3.07994089e-04,
-3.22560277e-04, -3.36962062e-04, -3.51125513e-04, -3.64869909e-04,
-3.78187035e-04, -3.91076892e-04, -4.03539479e-04, -4.15574797e-04,
-4.27182846e-04, -4.38363624e-04, -4.49117134e-04, -4.59443373e-04,
-4.69342344e-04, -4.78814044e-04, -4.87858475e-04, -4.96475637e-04,
-5.04665529e-04, -5.12428152e-04, -5.19763505e-04, -5.26671588e-04,
-5.33152402e-04, -5.39205947e-04, -5.44832222e-04, -5.50022519e-04,
-5.54724596e-04, -5.58929744e-04, -5.62637963e-04, -5.65849253e-04,
-5.68563614e-04, -5.70781047e-04, -5.72501551e-04, -5.73725126e-04,
-5.74451772e-04, -5.74681489e-04, -5.74414278e-04, -5.73650138e-04,
-5.72389069e-04, -5.70631071e-04, -5.68376144e-04, -5.65624289e-04,
-5.62375505e-04, -5.58629792e-04, -5.54387150e-04, -5.49647580e-04,
-5.44417699e-04, -5.38843111e-04, -5.32983382e-04, -5.26838512e-04,
-5.20408500e-04, -5.13693348e-04, -5.06693054e-04, -4.99407618e-04,
-4.91837041e-04, -4.83981323e-04, -4.75840464e-04, -4.67414463e-04,
-4.58703321e-04, -4.49707038e-04, -4.40425613e-04, -4.30859047e-04,
...
4.30859047e-04, 4.40425613e-04, 4.49707038e-04, 4.58703321e-04,
4.67414463e-04, 4.75840464e-04, 4.83981323e-04, 4.91837041e-04,
4.99407618e-04, 5.06693054e-04, 5.13693348e-04, 5.20408500e-04,
5.26838512e-04, 5.32983382e-04, 5.38843111e-04, 5.44417699e-04,
5.49647580e-04, 5.54387150e-04, 5.58629792e-04, 5.62375505e-04,
5.65624289e-04, 5.68376144e-04, 5.70631071e-04, 5.72389069e-04,
5.73650138e-04, 5.74414278e-04, 5.74681489e-04, 5.74451772e-04,
5.73725126e-04, 5.72501551e-04, 5.70781047e-04, 5.68563614e-04,
5.65849253e-04, 5.62637963e-04, 5.58929744e-04, 5.54724596e-04,
5.50022519e-04, 5.44832222e-04, 5.39205947e-04, 5.33152402e-04,
5.26671588e-04, 5.19763505e-04, 5.12428152e-04, 5.04665529e-04,
4.96475637e-04, 4.87858475e-04, 4.78814044e-04, 4.69342344e-04,
4.59443373e-04, 4.49117134e-04, 4.38363624e-04, 4.27182846e-04,
4.15574797e-04, 4.03539479e-04, 3.91076892e-04, 3.78187035e-04,
3.64869909e-04, 3.51125513e-04, 3.36962062e-04, 3.22560277e-04,
3.07994089e-04, 2.93263499e-04, 2.78368506e-04, 2.63309110e-04,
2.48085311e-04, 2.32697109e-04, 2.17144505e-04, 2.01427498e-04,
1.85546089e-04, 1.69500276e-04, 1.53290061e-04, 1.36915443e-04,
1.20376422e-04, 1.03672999e-04, 8.68051725e-05, 6.97729435e-05,
5.25763117e-05, 3.52152773e-05, 1.76898400e-05, 1.25504872e-16])
Coordinates:
t float64 8B 0.05
* eta2 (eta2) float64 2kB 0.0 0.003922 0.007843 ... 0.9922 0.9961 1.0
component int64 8B 0
t_seconds float64 8B 1.668e-10
X (eta2) float64 2kB 1.34 1.34 1.34 1.339 ... 1.339 1.34 1.34 1.34
Y (eta2) float64 2kB -0.0 -0.0 -0.0 -0.0 ... -0.0 -0.0 -0.0 -0.0
Z (eta2) float64 2kB 0.0 0.008377 0.01675 ... -0.008377 -8.328e-17
Attributes:
run: dt=0.05, algo=LieTrotter, Nel=(6, 12, 2), p=(2, 2, 2)
run_name: representation_torus
Vector field <xarray.DataArray 'e_field' (eta2: 256)> Size: 2kB
array([ 6.27905524e-16, -8.84492000e-05, -1.76076386e-04, -2.62881559e-04,
-3.48864717e-04, -4.34025862e-04, -5.18364994e-04, -6.01882111e-04,
-6.84577215e-04, -7.66450305e-04, -8.47501381e-04, -9.27730443e-04,
-1.00713749e-03, -1.08572253e-03, -1.16348555e-03, -1.24042655e-03,
-1.31654555e-03, -1.39184253e-03, -1.46631749e-03, -1.53997045e-03,
-1.61280138e-03, -1.68481031e-03, -1.75562756e-03, -1.82434954e-03,
-1.89093518e-03, -1.95538446e-03, -2.01769740e-03, -2.07787399e-03,
-2.13591423e-03, -2.19181812e-03, -2.24558567e-03, -2.29721687e-03,
-2.34671172e-03, -2.39407022e-03, -2.43929238e-03, -2.48237819e-03,
-2.52332765e-03, -2.56214076e-03, -2.59881752e-03, -2.63335794e-03,
-2.66576201e-03, -2.69602973e-03, -2.72416111e-03, -2.75011260e-03,
-2.77362298e-03, -2.79464872e-03, -2.81318981e-03, -2.82924627e-03,
-2.84281807e-03, -2.85390524e-03, -2.86250775e-03, -2.86862563e-03,
-2.87225886e-03, -2.87340745e-03, -2.87207139e-03, -2.86825069e-03,
-2.86194534e-03, -2.85315536e-03, -2.84188072e-03, -2.82812145e-03,
-2.81187752e-03, -2.79314896e-03, -2.77193575e-03, -2.74823790e-03,
-2.72208849e-03, -2.69421556e-03, -2.66491691e-03, -2.63419256e-03,
-2.60204250e-03, -2.56846674e-03, -2.53346527e-03, -2.49703809e-03,
-2.45918521e-03, -2.41990662e-03, -2.37920232e-03, -2.33707232e-03,
-2.29351661e-03, -2.24853519e-03, -2.20212807e-03, -2.15429524e-03,
...
2.15429524e-03, 2.20212807e-03, 2.24853519e-03, 2.29351661e-03,
2.33707232e-03, 2.37920232e-03, 2.41990662e-03, 2.45918521e-03,
2.49703809e-03, 2.53346527e-03, 2.56846674e-03, 2.60204250e-03,
2.63419256e-03, 2.66491691e-03, 2.69421556e-03, 2.72208849e-03,
2.74823790e-03, 2.77193575e-03, 2.79314896e-03, 2.81187752e-03,
2.82812145e-03, 2.84188072e-03, 2.85315536e-03, 2.86194534e-03,
2.86825069e-03, 2.87207139e-03, 2.87340745e-03, 2.87225886e-03,
2.86862563e-03, 2.86250775e-03, 2.85390524e-03, 2.84281807e-03,
2.82924627e-03, 2.81318981e-03, 2.79464872e-03, 2.77362298e-03,
2.75011260e-03, 2.72416111e-03, 2.69602973e-03, 2.66576201e-03,
2.63335794e-03, 2.59881752e-03, 2.56214076e-03, 2.52332765e-03,
2.48237819e-03, 2.43929238e-03, 2.39407022e-03, 2.34671172e-03,
2.29721687e-03, 2.24558567e-03, 2.19181812e-03, 2.13591423e-03,
2.07787399e-03, 2.01769740e-03, 1.95538446e-03, 1.89093518e-03,
1.82434954e-03, 1.75562756e-03, 1.68481031e-03, 1.61280138e-03,
1.53997045e-03, 1.46631749e-03, 1.39184253e-03, 1.31654555e-03,
1.24042655e-03, 1.16348555e-03, 1.08572253e-03, 1.00713749e-03,
9.27730443e-04, 8.47501381e-04, 7.66450305e-04, 6.84577215e-04,
6.01882111e-04, 5.18364994e-04, 4.34025862e-04, 3.48864717e-04,
2.62881559e-04, 1.76076386e-04, 8.84492000e-05, 6.27524359e-16])
Coordinates:
t float64 8B 0.05
* eta2 (eta2) float64 2kB 0.0 0.003922 0.007843 ... 0.9922 0.9961 1.0
component int64 8B 0
t_seconds float64 8B 1.668e-10
X (eta2) float64 2kB 1.34 1.34 1.34 1.339 ... 1.339 1.34 1.34 1.34
Y (eta2) float64 2kB -0.0 -0.0 -0.0 -0.0 ... -0.0 -0.0 -0.0 -0.0
Z (eta2) float64 2kB 0.0 0.008377 0.01675 ... -0.008377 -8.328e-17
Attributes:
run: dt=0.05, algo=LieTrotter, Nel=(6, 12, 2), p=(2, 2, 2)
run_name: representation_torus
Apply the workflow to another run#
For an already completed simulation, possibly in a separate process without MPI, open its output folder:
import struphy
out = struphy.Output("/path/to/sim_1").pproc(physical=True)
out.domain, out.model.units # reconstructed directly from saved metadata
Use out.scalars, out.fields, out.distributions, out.orbits, and out.densities. Attribute access is the normal interactive API; the corresponding *_catalog mappings are intended for generic loops and tooling.