Symplectic integration in a user-defined magnetostatic vector potential#964
Symplectic integration in a user-defined magnetostatic vector potential#964cemitch99 wants to merge 72 commits into
Conversation
for more information, see https://pre-commit.ci
for more information, see https://pre-commit.ci
for more information, see https://pre-commit.ci
for more information, see https://pre-commit.ci
9a1e4af to
fa61eba
Compare
…linear_element # Conflicts: # .github/workflows/cuda.yml # examples/CMakeLists.txt
Merging this PR will improve performance by 5.58%
|
| Benchmark | BASE |
HEAD |
Efficiency | |
|---|---|---|---|---|
| ⚡ | test_ExactDrift[spin] |
645.1 µs | 611 µs | +5.58% |
Tip
Curious why this is faster? Comment @codspeedbot explain why this is faster on this PR, or directly use the CodSpeed MCP with your agent.
Comparing cemitch99:add_user_nonlinear_element (0f461cf) with development (69c4e3e)
| * local frame. Reversing only negates the path-length integration | ||
| * direction; the user is responsible for the field's z-symmetry. | ||
| */ | ||
| void reverse () { Thick::reverse(); } |
There was a problem hiding this comment.
To check: complete or is there more to do? (needs a test entry in tests/python/test_reversibility_elements.py)
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
|
@cemitch99 I updated this PR and test pass now :) A few smaller inline comments remain. |
| // Iteration loop (implicit step) | ||
| for (int j = 0; j < max_iterations; j++) { | ||
|
|
There was a problem hiding this comment.
Codex:
- High: src/particles/integrators/Integrators.H:642 silently leaves particle
unchanged if the fixed-point/Newton loop never reaches tolerance. This can
turn a hard integration failure into a no-op push for that particle. I’d
abort/report non-convergence, and make the tolerance precision-aware.
There was a problem hiding this comment.
Could be tricky within a GPU kernel. We could make particles invalid with signaling NaNs or whip up a detailed error machinery.
| .def("to_dict", | ||
| [](MagnetostaticVectorPotential const & vp) { | ||
| return element_dict( | ||
| vp, | ||
| std::make_pair("unit", vp.m_unit) | ||
| ); | ||
| } | ||
| ) |
There was a problem hiding this comment.
Codex:
- Medium: src/python/elements.cpp:2886 serializes only unit for
MagnetostaticVectorPotential. ax, ay, derivative formulas, int_order, and
mapsteps are dropped, so KnownElementsList.to_dicts()/from_dicts()
reconstructs a different zero-field/default element. The new element should
also be added to tests/python/test_element_serialization.py.
| .. py:property:: k | ||
|
|
||
| quadrupole strength in 1/m^2 (or T/m) |
There was a problem hiding this comment.
Codex:
- Low: docs/source/usage/python.rst:1807 documents a nonexistent k property
and describes unit as a quad strength unit. Also mapsteps is documented as
default 5 in docs, while the C++/Python binding default is 10.
| void | ||
| finalize () | ||
| { | ||
| DynamicData::registry().erase(m_id); | ||
| } |
There was a problem hiding this comment.
Codex:
- Medium: src/elements/MagnetostaticVectorPotential.H:409 erases the parser
registry entry during finalize(). Since elements are copied by value into
the lattice, a user-held Python/C++ element object can still exist with the
same m_id after sim.finalize(), but its parser data is gone. Reusing that
element in a later simulation will fail in DynamicData::get(m_id). This also
contradicts the comment at src/elements/MagnetostaticVectorPotential.H:44
saying host data persists across simulations. A robust fix would keep the
original formula strings as durable host data and rebuild AMReX parsers
after a new initialize cycle.
There was a problem hiding this comment.
Yes this pattern is wrong, needs to be a NoFinalize class. Will update.
| } else { | ||
|
|
||
| // initialize numerical integration parameters (int_order = 6) |
There was a problem hiding this comment.
Codex:
- Medium: src/particles/integrators/Integrators.H:568 treats every int_order
other than 2 or 4 as sixth order. Because the constructor, input parser, and
Python setter accept arbitrary integers, int_order=3 or a typo silently
changes the integration method. This should be else if (int_order == 6) plus
an error, or validation in the element constructor/parser/setter.
There was a problem hiding this comment.
I agree, and the setter/constructor should abort on orders that are not 2,4 or 6
|
|
||
| .. literalinclude:: analysis_fodo_vector_potential.py | ||
| :language: python3 | ||
| :caption: You can copy this file from ``examples/vector_potential/analysis_vector_potential.py``. |
There was a problem hiding this comment.
Codex:
- Low: examples/vector_potential/README.rst:45 has a stale caption pointing to
analysis_vector_potential.py, but the included file is
analysis_fodo_vector_potential.py.
Co-authored-by: Chad Mitchell <46825199+cemitch99@users.noreply.github.com>
| element.map(tau2,particle2,particle1,zeval2,zeval1); | ||
| element.map(tau1,particle1,particle2,zeval1,zeval2); | ||
|
|
||
| } else { |
There was a problem hiding this comment.
| } else { | |
| } else if (int_order==6) { |
| element.map(tau3,particle1,particle2,zeval1,zeval2); | ||
| element.map(tau2,particle2,particle1,zeval2,zeval1); | ||
| element.map(tau1,particle1,particle2,zeval1,zeval2); | ||
|
|
There was a problem hiding this comment.
| } else { | |
| // throw an error here |
This PR adds a new element that allows the user to track through a region with a specified magnetostatic vector potential. Symplectic integration is performed using the exact form of the nonlinear relativistic Hamiltonian. We use the semiexplicit integrator appearing in:
B. Jayawardana and T. Ohsawa, "Semiexplicit symplectic integrators for non-separable Hamiltonian systems," Math. Comput. 92, pp. 251-281 (2022),
https://doi.org/10.1090/mcom/3778 (arXiv)
To do: