helicon.config

Configuration parsing, validation, and WarpX input file generation.

Configuration Classes

helicon.config.parser.SimConfig

Bases: BaseModel

Top-level simulation configuration.

Functions

from_preset(name) classmethod

Load a built-in preset configuration by name.

Available presets: sunbird, dfd, ppr.

Source code in src/helicon/config/parser.py

@classmethod
def from_preset(cls, name: str) -> SimConfig:
    """Load a built-in preset configuration by name.

    Available presets: ``sunbird``, ``dfd``, ``ppr``.
    """
    preset_dir = importlib.resources.files("helicon.config.presets")
    preset_file = preset_dir / f"{name}.yaml"
    if not preset_file.is_file():
        available = [
            p.name.removesuffix(".yaml")
            for p in preset_dir.iterdir()
            if p.name.endswith(".yaml")
        ]
        msg = f"Unknown preset {name!r}. Available: {available}"
        raise ValueError(msg)
    text = preset_file.read_text()
    data = yaml.safe_load(text)
    return cls.model_validate(data)

from_yaml(path) classmethod

Load configuration from a YAML file.

Source code in src/helicon/config/parser.py

@classmethod
def from_yaml(cls, path: str | Path) -> SimConfig:
    """Load configuration from a YAML file."""
    path = Path(path)
    with path.open() as f:
        data = yaml.safe_load(f)
    return cls.model_validate(data)

to_yaml(path)

Write this configuration to a YAML file.

Source code in src/helicon/config/parser.py

def to_yaml(self, path: str | Path) -> None:
    """Write this configuration to a YAML file."""
    path = Path(path)
    data = self.model_dump(mode="python")
    with path.open("w") as f:
        yaml.dump(data, f, default_flow_style=False, sort_keys=False)

helicon.config.parser.NozzleConfig

Bases: BaseModel

Nozzle geometry definition.

helicon.config.parser.CoilConfig

Bases: BaseModel

A single circular current loop.

helicon.config.parser.DomainConfig

Bases: BaseModel

Simulation domain extent.

helicon.config.parser.ResolutionConfig

Bases: BaseModel

Grid resolution.

helicon.config.parser.PlasmaSourceConfig

Bases: BaseModel

Upstream plasma injection boundary condition.

helicon.config.parser.DiagnosticsConfig

Bases: BaseModel

Diagnostic output configuration.

Validation

helicon.config.validators.validate_config(config)

Run all physics sanity checks on a simulation configuration.

Returns a ValidationResult with any warnings or errors found.

Source code in src/helicon/config/validators.py

def validate_config(config: SimConfig) -> ValidationResult:
    """Run all physics sanity checks on a simulation configuration.

    Returns a ``ValidationResult`` with any warnings or errors found.
    """
    warns: list[str] = []
    errors: list[str] = []

    # --- Domain checks ---
    domain = config.nozzle.domain
    domain_length = domain.z_max - domain.z_min
    if domain_length < 0.1:
        warns.append(
            f"Domain axial length ({domain_length:.3f} m) is very short. "
            "Detachment may not be captured."
        )

    # --- Coil checks ---
    for i, coil in enumerate(config.nozzle.coils):
        if coil.r > domain.r_max:
            errors.append(
                f"Coil {i} radius ({coil.r:.3f} m) exceeds domain r_max "
                f"({domain.r_max:.3f} m)."
            )
        if not (domain.z_min <= coil.z <= domain.z_max):
            warns.append(
                f"Coil {i} at z={coil.z:.3f} m is outside the domain "
                f"[{domain.z_min}, {domain.z_max}]."
            )

    # --- Resolution checks ---
    res = config.nozzle.resolution
    dr = domain.r_max / res.nr
    min_coil_r = min(c.r for c in config.nozzle.coils)
    if min_coil_r / dr < 5:
        warns.append(
            f"Radial resolution may be insufficient: smallest coil radius "
            f"({min_coil_r:.3f} m) is only {min_coil_r / dr:.1f} cells wide."
        )

    # --- Plasma beta estimate at throat ---
    # beta = n * k_B * (T_i + T_e) / (B^2 / 2 mu_0)
    # Rough B at throat: use on-axis field from strongest coil
    strongest = max(config.nozzle.coils, key=lambda c: abs(c.I))
    B_throat_approx = MU_0 * abs(strongest.I) / (2.0 * strongest.r)

    eV_to_J = 1.602176634e-19
    n0 = config.plasma.n0
    T_total_J = (config.plasma.T_i_eV + config.plasma.T_e_eV) * eV_to_J
    p_plasma = n0 * T_total_J
    p_mag = B_throat_approx**2 / (2.0 * MU_0)

    if p_mag > 0:
        beta = p_plasma / p_mag
        if beta > 1.0:
            warns.append(
                f"Estimated plasma beta at throat ≈ {beta:.2f} > 1. "
                "Magnetic confinement may be insufficient."
            )

    # --- Mass ratio warning ---
    if config.plasma.mass_ratio is not None and config.plasma.mass_ratio < 100:
        warns.append(
            f"Mass ratio {config.plasma.mass_ratio} is very low. "
            "Results will be qualitative only."
        )

    # --- Timestep check ---
    if config.timesteps < 1000:
        warns.append(f"Only {config.timesteps} timesteps — may not reach steady state.")

    # Emit warnings to Python warning system as well
    for w in warns:
        warnings.warn(w, UserWarning, stacklevel=2)

    return ValidationResult(
        passed=len(errors) == 0,
        warnings=warns,
        errors=errors,
    )

WarpX Input Generation

helicon.config.warpx_generator.generate_warpx_input(config)

Generate a WarpX input script string from a SimConfig.

Returns the input file content as a string. This can be written to disk and passed to WarpX via warpx.inputs.

Source code in src/helicon/config/warpx_generator.py

def generate_warpx_input(config: SimConfig) -> str:
    """Generate a WarpX input script string from a SimConfig.

    Returns the input file content as a string. This can be written
    to disk and passed to WarpX via ``warpx.inputs``.
    """
    domain = config.nozzle.domain
    res = config.nozzle.resolution
    cfg_hash = _config_hash(config)
    has_electrons = "e-" in config.plasma.species

    # Random seed
    seed = config.random_seed
    if seed is None:
        seed = int(cfg_hash, 16) % (2**31)

    is_3d = getattr(res, "geometry", "2d_rz") == "3d"

    lines = [
        f"# WarpX input file generated by Helicon v{_helicon_version()}",
        f"# Config hash: {cfg_hash}",
        "#",
        f"# Geometry: {'3D Cartesian' if is_3d else '2D RZ axisymmetric'}",
        f"# Nozzle type: {config.nozzle.type}",
        "",
        "################",
        "# GEOMETRY",
        "################",
    ]

    if is_3d:
        np_phi = getattr(res, "np_phi", 32)
        r_max = domain.r_max
        lines.extend(
            [
                "geometry.dims = 3",
                f"geometry.prob_lo = {domain.z_min} {-r_max} {-r_max}",
                f"geometry.prob_hi = {domain.z_max} {r_max} {r_max}",
                "",
                f"amr.n_cell = {res.nz} {res.nr} {np_phi}",
                "amr.max_level = 0",
                "",
                "# NOTE: 3D geometry is expensive on Apple Silicon CPU.",
                "# Use cloud HPC offload (helicon v1.3) for production 3D runs.",
                "",
                "################",
                "# BOUNDARY (3D)",
                "################",
                "boundary.field_lo = none none pec",
                "boundary.field_hi = pml pml pml",
                "boundary.particle_lo = absorbing absorbing absorbing",
                "boundary.particle_hi = absorbing absorbing absorbing",
            ]
        )
    else:
        lines.extend(
            [
                "geometry.dims = RZ",
                f"geometry.prob_lo = {domain.z_min} 0.0",
                f"geometry.prob_hi = {domain.z_max} {domain.r_max}",
                "",
                f"amr.n_cell = {res.nz} {res.nr}",
                "amr.max_level = 0",
                "",
                "################",
                "# BOUNDARY",
                "################",
                "boundary.field_lo = none pec",
                "boundary.field_hi = pml pml",
                "boundary.particle_lo = absorbing reflecting",
                "boundary.particle_hi = absorbing absorbing",
            ]
        )

    lines.extend(
        [
            "",
            "################",
            "# TIMESTEP",
            "################",
            f"max_step = {config.timesteps}",
            f"warpx.cfl = {config.dt_multiplier}",
            "",
            "################",
            "# PARTICLES",
            "################",
        ]
    )

    species_names = []
    for sp in config.plasma.species:
        safe_name = sp.replace("+", "_plus").replace("-", "_minus")
        species_names.append(safe_name)

    lines.append(f"particles.species_names = {' '.join(species_names)}")
    lines.append("")

    eV_to_J = 1.602176634e-19
    m_e_physical = SPECIES_MASS["e-"]

    # Reduced mass ratio: override electron mass
    reduced_mass_ratio = config.plasma.mass_ratio is not None and has_electrons
    if reduced_mass_ratio:
        mr = config.plasma.mass_ratio
        # Quantitative guidance (spec §8 v0.3: reduced mass ratio guidance documentation)
        # Physical D/e ratio = 3670.  Error bounds from PIC benchmarking literature:
        #   mr < 100  → qualitative only; η_d errors typically >30 %
        #   100 ≤ mr < 500  → qualitative trends correct; ~20–30 % η_d error
        #   500 ≤ mr < 1836 → semi-quantitative; ~5–15 % η_d error
        #   mr ≥ 1836       → near-physical; errors typically <5 % for D/e
        if mr < 100:
            _mr_guidance = (
                "Qualitative trends only; quantitative errors >30 % in η_d. "
                "Do not cite detachment numbers from this run."
            )
        elif mr < 500:
            _mr_guidance = (
                f"~20–30 % error in η_d vs full ratio ({3670:.0f}). "
                "Suitable for fast parameter-space exploration only."
            )
        elif mr < 1836:
            _mr_guidance = (
                f"~5–15 % error in η_d vs full ratio ({3670:.0f}). "
                "Semi-quantitative; use full ratio for validation cases."
            )
        else:
            _mr_guidance = f"Near-physical ratio; η_d errors typically <5 % vs {3670:.0f}."
        log.warning(
            "Reduced mass ratio m_i/m_e = %.1f requested. %s "
            "Output will be flagged mass_ratio_reduced=true.",
            mr,
            _mr_guidance,
        )
        # Use the lightest ion mass to compute effective electron mass
        ion_species = [s for s in config.plasma.species if s != "e-"]
        _default_mass = 1.67e-27
        if ion_species:
            m_ion_ref = min(SPECIES_MASS.get(s, _default_mass) for s in ion_species)
        else:
            m_ion_ref = _default_mass
        m_e_effective = m_ion_ref / mr
    else:
        m_e_effective = m_e_physical

    # Fluid-hybrid electron mode
    fluid_electrons = config.plasma.electron_model == "fluid" and has_electrons
    if fluid_electrons:
        lines.extend(
            [
                "################",
                "# HYBRID ELECTRON MODEL",
                "################",
                "# Fluid (hybrid) electron model: faster than fully kinetic,",
                "# less accurate for electron detachment physics.",
                "warpx.do_hybrid_model = 1",
                f"hybrid_model.Te_eV = {config.plasma.T_e_eV:.4f}",
                f"hybrid_model.n0_m3 = {config.plasma.n0:.6e}",
                "hybrid_model.plasma_resistivity = 0.0",
                "",
            ]
        )

    for sp, safe_name in zip(config.plasma.species, species_names, strict=True):
        # Skip electron particle species when using fluid model
        if sp == "e-" and fluid_electrons:
            continue

        mass = SPECIES_MASS.get(sp)
        charge = SPECIES_CHARGE.get(sp)
        if mass is None or charge is None:
            lines.append(f"# WARNING: Unknown species {sp}, using placeholder values")
            mass = mass or 1.67262192595e-27
            charge = charge or 1.602176634e-19

        # Override electron mass for reduced mass ratio
        if sp == "e-" and reduced_mass_ratio:
            mass = m_e_effective

        T_eV = config.plasma.T_e_eV if sp == "e-" else config.plasma.T_i_eV
        T_J = T_eV * eV_to_J
        v_th = math.sqrt(T_J / mass)

        lines.extend(
            [
                (
                    f"{safe_name}.species_type = proton"
                    if "+" in sp
                    else f"{safe_name}.species_type = electron"
                ),
                f"{safe_name}.mass = {mass:.10e}",
                f"{safe_name}.charge = {charge:.10e}",
                f"{safe_name}.injection_style = NUniformPerCell",
                f"{safe_name}.num_particles_per_cell = 64",
                f"{safe_name}.profile = constant",
                f"{safe_name}.density = {config.plasma.n0:.6e}",
                f"{safe_name}.momentum_distribution_type = gaussian",
                f"{safe_name}.uz_m = {config.plasma.v_injection_ms:.6e}",
                f"{safe_name}.ux_th = {v_th:.6e}",
                f"{safe_name}.uy_th = {v_th:.6e}",
                f"{safe_name}.uz_th = {v_th:.6e}",
                "",
            ]
        )

    if reduced_mass_ratio:
        lines.append(
            f"# METADATA: mass_ratio_reduced = true (m_i/m_e = {config.plasma.mass_ratio})"
        )
    if fluid_electrons:
        lines.append("# METADATA: electron_model = fluid")

    # Monte Carlo neutrals (spec §2.1)
    neutrals = config.plasma.neutrals
    if neutrals is not None:
        neutral_name = neutrals.species.replace("+", "_plus").replace("-", "_minus")
        lines.extend(
            [
                "",
                "################",
                "# MONTE CARLO NEUTRALS",
                "################",
                f"# Background {neutrals.species} neutral gas: Monte Carlo Collision (MCC)",
                f"particles.species_names += {neutral_name}",
                f"{neutral_name}.species_type = neutral",
                f"{neutral_name}.mass = "
                f"{SPECIES_MASS.get(neutrals.species + '+', 3.34e-27):.10e}",
                f"{neutral_name}.charge = 0.0",
                f"{neutral_name}.injection_style = NUniformPerCell",
                f"{neutral_name}.num_particles_per_cell = 16",
                f"{neutral_name}.profile = constant",
                f"{neutral_name}.density = {neutrals.n_neutral_m3:.6e}",
                f"{neutral_name}.momentum_distribution_type = gaussian",
                f"{neutral_name}.ux_th = {_neutral_vth(neutrals, eV_to_J):.6e}",
                f"{neutral_name}.uy_th = {_neutral_vth(neutrals, eV_to_J):.6e}",
                f"{neutral_name}.uz_th = {_neutral_vth(neutrals, eV_to_J):.6e}",
                "",
                "# Charge-exchange collisions",
                "mcc.species = "
                + " ".join(
                    s.replace("+", "_plus").replace("-", "_minus")
                    for s in config.plasma.species
                    if "+" in s
                ),
                f"mcc.neutral_species = {neutral_name}",
                f"mcc.CX_cross_section = {neutrals.cx_cross_section_m2:.6e}",
            ]
        )
        if neutrals.ionization_cross_section_m2 is not None:
            lines.extend(
                [
                    f"mcc.ionization_cross_section = "
                    f"{neutrals.ionization_cross_section_m2:.6e}",
                    "mcc.do_ionization = 1",
                ]
            )

    lines.extend(
        [
            "################",
            "# EXTERNAL B-FIELD",
            "################",
            "# Applied B-field is loaded from HDF5 file pre-computed by Helicon",
            "# Use: warpx.read_fields_from_path = <bfield.h5>",
            "",
            "################",
            "# DIAGNOSTICS",
            "################",
        ]
    )

    # Build diags_names list once (avoids duplicate key)
    diag_names = ["diag1"]
    if config.diagnostics.mode == "analysis":
        diag_names.append("diag_particles")

    lines.append(f"diagnostics.diags_names = {' '.join(diag_names)}")
    lines.extend(
        [
            f"diag1.intervals = {config.diagnostics.field_dump_interval}",
            "diag1.diag_type = Full",
            "diag1.format = openpmd",
            "diag1.openpmd_backend = h5",
            "diag1.fields_to_plot = Br Bz Er Ez jr jz rho",
            "",
        ]
    )

    if config.diagnostics.mode == "analysis":
        lines.extend(
            [
                f"diag_particles.intervals = {config.diagnostics.particle_dump_interval}",
                "diag_particles.diag_type = Full",
                "diag_particles.format = openpmd",
                "diag_particles.openpmd_backend = h5",
                f"diag_particles.species = {' '.join(species_names)}",
                "",
            ]
        )

    lines.extend(
        [
            "################",
            "# REPRODUCIBILITY",
            "################",
            f"warpx.random_seed = {seed}",
            "warpx.verbose = 1",
        ]
    )

    return "\n".join(lines) + "\n"