emout

emout

A Python library for parsing and visualizing output files generated by EMSES simulations.

• Documentation: https://nkzono99.github.io/emout/

Installation

pip install emout

Example

• Visualization of simulation results for lunar surface charging

Overview

When you run EMSES simulations, the results (e.g., potentials, densities, currents) are output in .h5 files, and a parameter file (plasma.inp) contains the simulation settings. emout helps you:

• emout

  • Installation

  • Example

  • Overview

  • Usage

    • Loading Data

    • Retrieving the Parameter File (plasma.inp)

    • Plotting Data

    • Working with Units

    • Handling Appended Simulation Outputs

    • Data Masking

    • Creating Animations

    • Solving Poisson’s Equation (Experimental)

    • Backtrace Usage Examples (Experimental)

      • Install vdist-solver-fortran

      • Running backtrace on HPC computational nodes with Dask (>= Python3.10)

      • Backtrace using the particles from a previous probability calculation, then plot $x$–$z$ trajectories with probability as transparency

Below, you will find usage examples that assume the following directory structure:

.
└── output_dir
    ├── plasma.inp
    ├── phisp00_0000.h5
    ├── nd1p00_0000.h5
    ├── nd2p00_0000.h5
    ├── j1x00_0000.h5
    ├── j1y00_0000.h5
    ...
    └── bz00_0000.h5

---

Usage

Loading Data

import emout

# Initialize Emout with the path to the output directory
data = emout.Emout('output_dir')

# Access arrays by their variable names (derived from EMSES filename)
data.phisp   # Data from phisp00_0000.h5
len(data.phisp)   # Number of time steps
data.phisp[0].shape

data.j1x
data.bz
data.j1xy  # Vector data from "j1x00_0000.h5" + "j1y00_0000.h5"

# Data created by "relocating" ex00_0000.h5
data.rex

# Access data as a pandas DataFrame
data.icur
data.pbody

---

Retrieving the Parameter File (plasma.inp)

# The namelist is parsed into a dictionary-like structure
data.inp
data.inp['tmgrid']['nx']  # Access via group name and parameter name
data.inp['nx']            # Group name can be omitted if not ambiguous
data.inp.tmgrid.nx        # Access like object attributes
data.inp.nx               # Still valid

---

Plotting Data

x, y, z = 32, 32, 100

# Basic 2D plot (xy-plane at z=100 for the last timestep)
data.phisp[-1, z, :, :].plot()

# Line plot along z-axis at x=32, y=32
data.phisp[-1, :, y, x].plot()

# Plot using SI units
data.phisp[-1, z, :, :].plot(use_si=True)  # Default is True

# Show or save plot
data.phisp[-1, z, :, :].plot(show=True)
data.phisp[-1, z, :, :].plot(savefilename='phisp.png')

# Plot vector field as a streamline
data.j1xy[-1, z, :, :].plot()

---

Working with Units

Note EMSES → SI conversion is supported only when the first line of plasma.inp includes something like:

!!key dx=[0.5],to_c=[10000.0]

where dx is the grid spacing [m], and to_c is the (normalized) speed of light used internally by EMSES.

# Converting between physical (SI) and EMSES units
data.unit.v.trans(1)    # Convert 1 m/s to EMSES velocity unit
data.unit.v.reverse(1)  # Convert 1 EMSES velocity unit to m/s

# Access converted data in SI units
phisp_volt = data.phisp[-1, :, :, :].val_si       # Potential in volts [V]
j1z_A_per_m2 = data.j1z[-1, :, :, :].val_si       # Current density [A/m^2]
nd1p_per_cc = data.nd1p[-1, :, :, :].val_si       # Number density [1/cm^3]

Unit Name List

B = Magnetic flux density [T]
C = Capacitance [F]
E = Electric field [V/m]
F = Force [N]
G = Conductance [S]
J = Current density [A/m^2]
L = Inductance [H]
N = Flux [/m^2s]
P = Power [W]
T = Temperature [K]
W = Energy [J]
a = Acceleration [m/s^2]
c = Light Speed [m/s]
e = Napiers constant []
e0 = FS-Permttivity [F/m]
eps = Permittivity  [F/m]
f = Frequency [Hz]
i = Current [A]
kB = Boltzmann constant [J/K]
length = Sim-to-Real length ratio [m]
m = Mass [kg]
m0 = FS-Permeablity [N/A^2]
mu = Permiability [H/m]
n = Number density [/m^3]
phi = Potential [V]
pi = Circular constant []
q = Charge [C]
q_m = Charge-to-mass ratio [C/kg]
qe = Elementary charge [C]
qe_me = Electron charge-to-mass ratio [C/kg]
rho = Charge density [C/m^3]
t = Time [s]
v = Velocity [m/s]
w = Energy density [J/m^3]

---

Handling Appended Simulation Outputs

Examples

If your simulation continues and creates new directories:

import emout

# Merge multiple output directories into one Emout object
data = emout.Emout('output_dir', append_directories=['output_dir_2', 'output_dir_3'])

# Same as above if 'ad="auto"' is specified (detects appended outputs automatically)
data = emout.Emout('output_dir', ad='auto')

---

Data Masking

Examples

# Mask values below the average
data.phisp[1].masked(lambda phi: phi < phi.mean())

# Equivalent manual approach
phi = data.phisp[1].copy()
phi[phi < phi.mean()] = float('nan')

---

Creating Animations

Examples

# Create a time-series animation along the first axis (time = 0)
x, y, z = 32, 32, 100
data.phisp[:, z, :, :].gifplot()

# Specify a different axis (default is axis=0)
data.phisp[:, z, :, :].gifplot(axis=0)

# Save animation as a GIF
data.phisp[:, z, :, :].gifplot(action='save', filename='phisp.gif')

# Display the animation inline in a Jupyter notebook
data.phisp[:, z, :, :].gifplot(action='to_html')

# Combining multiple frames for a single animation
updater0 = data.phisp[:, z, :, :].gifplot(action='frames', mode='cmap')
updater1 = data.phisp[:, z, :, :].build_frame_updater(mode='cont')
updater2 = data.nd1p[:, z, :, :].build_frame_updater(mode='cmap', vmin=1e-3, vmax=20, norm='log')
updater3 = data.nd2p[:, z, :, :].build_frame_updater(mode='cmap', vmin=1e-3, vmax=20, norm='log')
updater4 = data.j2xy[:, z, :, :].build_frame_updater(mode='stream')

layout = [
    [
        [updater0, updater1],
        [updater2],
        [updater3, updater4]
    ]
]
animator = updater0.to_animator(layout=layout)
animator.plot(action='to_html')  # or 'save', 'show', etc.

---

Solving Poisson’s Equation (Experimental)

You can solve Poisson’s equation from 3D charge distributions using emout.poisson (depends on scipy):

Examples

import numpy as np
import scipy.constants as cn
from emout import Emout
from emout.utils import poisson

data = Emout('output_dir')
dx = data.inp.dx  # [m] Grid spacing
rho = data.rho[-1].val_si  # [C/m^3] Charge distribution
btypes = ["pdn"[i] for i in data.inp.mtd_vbnd]  # Boundary conditions

# Solve Poisson’s equation for potential
phisp = poisson(rho, dx=dx, btypes=btypes, epsilon_0=cn.epsilon_0)

# Compare with EMSES potential
np.allclose(phisp, data.phisp[-1])  # Should be True (within numerical tolerance)

---

Backtrace Usage Examples (Experimental)

Examples

Install vdist-solver-fortran

pip install git+https://github.com/Nkzono99/vdist-solver-fortran.git

Below are three example workflows demonstrating how to use the data.backtrace interface. All examples assume you have already created an Emout object named data.

with Dask

Running backtrace on HPC computational nodes with Dask (>= Python3.10)

If you’ve set up a Dask cluster via emout.distributed, all of the data.backtrace calls below will actually run on your computational nodes instead of your login node.

from emout.distributed import start_cluster, stop_cluster
import emout

# ① Dask クラスタを起動（SLURM ジョブを一時的に作成して Worker を常駐させる）
client = start_cluster(
    partition="gr20001a",   # 使用するキュー
    processes=1,            # プロセス数
    cores=112,              # コア数
    memory="60G",           # メモリ
    walltime="03:00:00",    # 最大実行時間
    scheduler_ip=None,  # ログインノード上の Scheduler IP (e.g. "10.10.64.1", Noneで自動検索)
    scheduler_port=32332,       # Scheduler ポート
)

# ② 通常の data.backtrace API を呼び出すだけで、
#    図のようにバックトレース関数群が計算ノード上で実行されます
data = emout.Emout("output_dir")
result = data.backtrace.get_probabilities(
    128, 128, 200,
    (-data.inp.path[0]*3, data.inp.path[0]*3, 500),
    1,
    (-data.inp.path[0]*4, data.inp.path[0]*3, 500),
    ispec=0,
    istep=-1,
    dt=data.inp.dt,
    max_step=100000,
    n_threads=112,
)
result.vxvz.plot()

# ③ 終了後はクライアントを閉じて Scheduler を停止
stop_cluster()

Backtrace using the particles from a previous probability calculation, then plot $x$–$z$ trajectories with probability as transparency

We take the particles array produced internally by get_probabilities(...), run backtraces on each of those particles, and then plot the $x–z$ projections of all backtraced trajectories.

We normalize each trajectory’s probability to the maximum probability across all phase‐grid cells, and pass that normalized array to alpha, so that high‐probability trajectories appear more opaque and low‐probability trajectories more transparent.

import matplotlib.pyplot as plt
import numpy as np
import emout

data = emout.Emout()

ispec = 0 # e.g. 0: electron, 1: ion, 2: photoelectron

# 1) Create ProbabilityResult
probability_result = data.backtrace.get_probabilities(
    128,
    128,
    60,
    (-data.inp.path[0] * 3, data.inp.path[0] * 3, 10),
    0,
    (-data.inp.path[0] * 3, 0, 10),
    ispec=ispec,
#   dt=data.inp.dt, # Set dt=-data.inp.dt if you want to forward-trace
)

# 2) Plot probability distribution
probability_result.vxvz.plot()

# 3) Extract the `particles` array and their associated `probabilities`
particles = probability_result.particles        # Sequence of Particle objects
prob_flat  = probability_result.probabilities    # 2D array of shape (nvz, nvx)

# 4) Flatten the 2D probability grid back into the 1D array matching `particles` order
prob_1d = prob_flat.ravel()

# 5) Normalize probabilities to [0,1] by dividing by the global maximum
alpha_values = np.nan_to_num(prob_1d / prob_1d.max())

# 6) Compute backtraces for all particles
backtrace_result = data.backtrace.get_backtraces_from_particles(
    particles,
    ispec=ispec,
#   dt=data.inp.dt, # Set dt=-data.inp.dt if you want to forward-trace
)

# 7) Plot x vs z for every trajectory, using the normalized probabilities as alpha
ax = backtrace_result.xz.plot(color="black", alpha=alpha_values)
ax.set_title("Backtrace Trajectories (x vs z) with Probability Transparency")
plt.show()

Notes on the above examples:

• get_backtraces(positions, velocities) returns a MultiBacktraceResult whose xy property is a MultiXYData object. You can sample, reorder, or subset the trajectories and then call .plot() on .xy, .vxvy, .xz, etc.

• get_probabilities(...) returns a ProbabilityResult whose .vxvz, .xy, .xz, etc. are all HeatmapData objects. Calling .plot() on any of these displays a 2D probability heatmap for the chosen pair of axes.

• probability_result.particles is the list of Particle objects used internally to compute the 6D probability grid. We pass that list to get_backtraces_from_particles(...) to compute backtraced trajectories for exactly those same particles. Normalizing their probabilities to [0,1] and passing that array into alpha makes high‐probability trajectories draw more opaquely.

These patterns demonstrate the flexibility of the data.backtrace facade for:

1. Direct backtracing from arbitrary $(\mathbf{r}, \mathbf{v})$ arrays,

2. Probability‐space calculations on a structured phase grid, and

3. Combining the two so that you can visualize backtraced trajectories with opacity weighted by their computed probabilities.

---

Feel free to explore the documentation for more details and advanced usage:

emout Documentation

Happy analyzing!

API Reference

• emout package
  • Subpackages
  • Module contents