Skip to content

161 add binning of higher order moments of f and delta f#162

Merged
spossann merged 27 commits intodevelfrom
161-add-binning-of-higher-order-moments-of-f-and-delta-f
Jan 28, 2026
Merged

161 add binning of higher order moments of f and delta f#162
spossann merged 27 commits intodevelfrom
161-add-binning-of-higher-order-moments-of-f-and-delta-f

Conversation

@P-Attapon
Copy link
Copy Markdown
Contributor

@P-Attapon P-Attapon commented Jan 16, 2026

Solves the following issue(s):
Add options to bin current density and energy tensor in BinningPlot instead of only particle density

Core changes:

  1. Implement weights calculation based on type of quantity to calculate.
  2. For saving data, the type of quantity is attached to the slice name.
    Example: slice = "e1_v1", output_quantity = "current_density_v1"
    => simdata.f["kinetic_ions"]["e1_v1_current_density_v1"]
    slice = "e1_v1", output_quantity = "particle"
    => simdata.f["kinetic_ions"]["e1_v1"]
  3. We make the signature and docstring of Particles base class appear in VScode, when an instance of Particle class is created:
    • new abstract method __post_init__(self) is called at the end of base __init__(self)
    • The handling of the default column numbers had to be changed
    • Move docstring from class to __init__

Model-specific changes:

None
Older parameter files can be compiled directly without any changes

Documentation changes:

  1. Create a new set of Literal BinningOutput. This Literal is used to determine the quantity to bin and the axis of velocity for weights in binning calculation. The three quantities are current_density, energy_tensor, and particle.

BinningOutput is used to determine quantity to be binned and respective velocity axis for weights calculation
struphy.pic.base.binning is edited. The weights are calculated based on 'output_quantity' defined in BinningOutput.
Edited to incoporate other binning plots besides particle density.
Binning different quntity in the same space is saved as a new key in simdata.f["kinetic_ions"].

Example:
binplot1 = BinningPlot(slice="e1_v1", ..., output_quantity="particle")
binplot2 = BinningPlot(slice="e1_v1", ..., output_quantity="current_density_v1")

simdata.f["kinetic_ions"].keys() => ["e1_v1", "e1_v1_current_density_v1"]
@P-Attapon P-Attapon requested a review from spossann January 16, 2026 14:45
@P-Attapon P-Attapon linked an issue Jan 16, 2026 that may be closed by this pull request
@P-Attapon P-Attapon marked this pull request as ready for review January 16, 2026 14:50
Copy link
Copy Markdown
Member

@spossann spossann left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This look pretty good! I have added some minor suggestions and the folllowing points:

  1. Please add unit test for the new binning outputs in pic/tests/test_binning. You can write new tests such as def test_binning_current_6D_full_f(mapping, show_plot=False): and adapt the following code accordingly:
# =========================================
    # ===== Test cosine in eta1 direction =====
    # =========================================
    # test weights
    amp_n = 0.1
    l_n = 2
    pert = perturbations.ModesCos(ls=(l_n,), amps=(amp_n,))
    maxwellian = Maxwellian3D(n=(1.0, pert))

    particles = Particles6D(
        loading_params=loading_params,
        boundary_params=boundary_params,
        domain=domain,
        background=maxwellian,
    )
    particles.draw_markers()
    particles.initialize_weights()

    e1_bins = xp.linspace(0.0, 1.0, 200, endpoint=True)
    de = e1_bins[1] - e1_bins[0]

    binned_res, r2 = particles.binning(
        [True, False, False, False, False, False],
        [e1_bins],
    )

    e1_plot = e1_bins[:-1] + de / 2

    ana_res = 1.0 + amp_n * xp.cos(2 * xp.pi * l_n * e1_plot)

    if show_plot:
        plt.plot(e1_plot, ana_res, label="Analytical result")
        plt.plot(e1_plot, binned_res, "r*", label="From binning")
        plt.title(r"Full-$f$: Cosine in $\eta_1$-direction")
        plt.xlabel(r"$\eta_1$")
        plt.ylabel(r"$f(\eta_1)$")
        plt.legend()
        plt.show()

    l2_error = xp.sqrt(xp.sum((ana_res - binned_res) ** 2)) / xp.sqrt(xp.sum((ana_res) ** 2))

    assert l2_error <= 0.02, f"Error between binned data and analytical result was {l2_error}"
  1. Please fix the unit test https://github.com/struphy-hub/struphy/actions/runs/21075207400/job/60615217672?pr=162#step:10:6877

  2. You have to run struphy format all to format the code correctly, such that the CI tests pass.

Comment thread src/struphy/io/options.py Outdated
Comment thread src/struphy/models/base.py Outdated
Comment thread src/struphy/pic/base.py Outdated
Comment thread src/struphy/pic/base.py Outdated
Comment thread src/struphy/pic/base.py Outdated
Comment thread src/struphy/pic/base.py Outdated
Comment thread src/struphy/pic/utilities.py Outdated
P-Attapon and others added 5 commits January 19, 2026 09:51
Co-authored-by: Stefan Possanner <86720346+spossann@users.noreply.github.com>
Rename literals of BinningQuantity to new convention.
if statements Particles.binning are cleand up
@P-Attapon
Copy link
Copy Markdown
Contributor Author

@spossann

Please take a look into the new binning for current

Also I could not do struphy format all:

struphy format all
usage: struphy [-h] [-v] [-s] [--fluid] [--kinetic] [--hybrid] [--toy] [--refresh-models] COMMAND ...
struphy: error: argument COMMAND: invalid choice: 'format' (choose from 'compile', 'params', 'profile', 'test')

Copy link
Copy Markdown
Member

@spossann spossann left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In order to be able to format you need to install as a developer pip install -e .[dev].

Comment thread src/struphy/pic/tests/test_binning.py Outdated
@P-Attapon P-Attapon marked this pull request as draft January 23, 2026 11:02
@P-Attapon P-Attapon marked this pull request as ready for review January 26, 2026 08:15
@P-Attapon
Copy link
Copy Markdown
Contributor Author

@spossann I've added all the test cases for the new binning we have made.

For the heat flux, the initial value does not depend on the perturbation in density, so the value is always zero. Let me know if I should change the perturbation to be in the velocity, in that case I am not sure how the analytical result should be calculated.

Best,

when instance of Particle class is created:

* new abstract method `__post_init__(self)` is called at the end of base
  __init__
* The handling of the default column numbers had to be changed
* Move docstring from class to __init__
spossann
spossann previously approved these changes Jan 26, 2026
Copy link
Copy Markdown
Member

@spossann spossann left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's fine now. I have added a small restructuring of the Particles base class for better argument info in VScode. We can adapt the test for the heat flux later, when actually needed.

Comment thread src/struphy/pic/base.py Outdated
Comment thread src/struphy/pic/tests/test_binning.py
@spossann
Copy link
Copy Markdown
Member

spossann commented Jan 26, 2026

@max-models our local formatting does not catch the ill-formatted test_binning.py. Do you have an idea what to do in such a situation? We ran struphy format all.

@max-models
Copy link
Copy Markdown
Member

@max-models our local formatting does not catch the ill-formatted test_binning.py. Do you have an idea what to do in such a situation? We ran struphy format all.

A quick fix will be to run isort ., but I'm not sure why there is a difference between ruff and isort here, seems like it has to do with extra newlines.

@spossann spossann merged commit ad45f3b into devel Jan 28, 2026
12 checks passed
@spossann spossann deleted the 161-add-binning-of-higher-order-moments-of-f-and-delta-f branch January 28, 2026 07:03
max-models pushed a commit that referenced this pull request Mar 20, 2026
**Solves the following issue(s):**
Add options to bin **current density** and **energy tensor** in
`BinningPlot` instead of only particle density

**Core changes:**

1. Implement weights calculation based on type of quantity to calculate.
2. For saving data, the type of quantity is attached to the slice name.
Example: slice = "e1_v1", output_quantity = "current_density_v1"
=> `simdata.f["kinetic_ions"]["e1_v1_current_density_v1"]`
slice = "e1_v1", output_quantity = "particle"
=>` simdata.f["kinetic_ions"]["e1_v1"]`
3. We make the signature and docstring of `Particles` base class appear
in VScode, when an instance of Particle class is created:
* new abstract method `__post_init__(self)` is called at the end of base
`__init__(self)`
    * The handling of the default column numbers had to be changed
    * Move docstring from class to `__init__`

**Model-specific changes:**

None
Older parameter files can be compiled directly without any changes

**Documentation changes:**

1. Create a new set of Literal `BinningOutput`. This Literal is used to
determine the quantity to bin and the axis of velocity for weights in
binning calculation. The three quantities are `current_density`,
`energy_tensor`, and `particle`.

---------

Co-authored-by: Stefan Possanner <86720346+spossann@users.noreply.github.com>
Co-authored-by: Stefan Possanner <stefan.possanner@ipp.mpg.de>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Add binning of higher-order moments of f and delta-f

3 participants