From 4d73a34820682895831804c56ffe0dbc066452ce Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Mon, 19 Oct 2026 17:18:31 +0200 Subject: [PATCH 1/2] DriftKineticElectrostaticAdiabatic: Poisson solve between the two guiding-center pushers Reorder the propagators from [gc_poisson, push_gc_bxe, push_gc_para] to [push_gc_bxe, gc_poisson, push_gc_para]. Both pushers move the guiding centers, so with Strang splitting the old order pushed para(dt) and the closing bxe(dt/2) with a potential solved before a drift. Now the Strang step is bxe(dt/2), P, para(dt), P, bxe(dt/2) at the same cost of two Poisson solves. Also solve the Poisson equation in update_scalar_quantities (once per saved step), so that en_phi, en_tot and the saved phi belong to the current marker positions. Co-Authored-By: Claude Opus 5.5 --- .../drift_kinetic_electrostatic_adiabatic.py | 37 +++++++++++++++---- 1 file changed, 29 insertions(+), 8 deletions(-) diff --git a/src/struphy/models/drift_kinetic_electrostatic_adiabatic.py b/src/struphy/models/drift_kinetic_electrostatic_adiabatic.py index 437ac7780..ff70059c2 100644 --- a/src/struphy/models/drift_kinetic_electrostatic_adiabatic.py +++ b/src/struphy/models/drift_kinetic_electrostatic_adiabatic.py @@ -89,8 +89,11 @@ def __init__( Z: int = 1, diagnostic: FEECVariable | None = None, ): - self.gc_poisson = GyrokineticPoissonSolve(rho=rho, epsilon=epsilon, Z=Z, diagnostic=diagnostic) + # both pushers move the guiding centers; with the Poisson solve between them, a Strang step + # bxe(dt/2), poisson, para(dt), poisson, bxe(dt/2) pushes with the potential of the current positions + # (the two bxe(dt/2) of consecutive steps act as one bxe(dt) with the potential frozen) self.push_gc_bxe = PushGuidingCenterBxEstar(phi=phi) + self.gc_poisson = GyrokineticPoissonSolve(rho=rho, epsilon=epsilon, Z=Z, diagnostic=diagnostic) self.push_gc_para = PushGuidingCenterParallel(phi=phi) ## abstract methods @@ -138,8 +141,8 @@ def __init__( ) # 4. assign variables to propagators - self.propagators.gc_poisson.variables.phi = self.em_fields.phi self.propagators.push_gc_bxe.variables.ions = self.kinetic_ions.var + self.propagators.gc_poisson.variables.phi = self.em_fields.phi self.propagators.push_gc_para.variables.ions = self.kinetic_ions.var # 5. define scalars to be tracked during simulation @@ -183,6 +186,18 @@ def post_allocate(self): average_matrix_teta = AverageOperator(derham, direction=1) self.average_matrix = average_matrix_z @ average_matrix_teta + def update_scalar_quantities(self): + """Solve the Poisson equation for the current marker positions, then update the scalars. + + A time step ends with a push (bxe with Strang, para with Lie-Trotter), so the potential would lag behind the + markers in ``en_phi``, ``en_tot`` and the saved data (which is extracted after this call). The extra solve, + once per saved step, does not change the dynamics: the next Poisson solve of the time stepping overwrites it. + + :meta private: + """ + self.propagators.gc_poisson(1.0) + super().update_scalar_quantities() + def _compute_en_phi(self): phi = self.em_fields.phi.spline.vector epsilon = self.kinetic_ions.equation_params.epsilon @@ -277,18 +292,24 @@ def doc_scalar_quantities(cls): def doc_discretization(cls): """Time integration is performed by the following propagators (in sequence): - 1. :class:`~struphy.propagators.implicit_diffusion.ImplicitDiffusion` - 2. :class:`~struphy.propagators.push_guiding_center_bx_estar.PushGuidingCenterBxEstar` + 1. :class:`~struphy.propagators.push_guiding_center_bx_estar.PushGuidingCenterBxEstar` + 2. :class:`~struphy.propagators.gyrokinetic_poisson_solve.GyrokineticPoissonSolve` 3. :class:`~struphy.propagators.push_guiding_center_parallel.PushGuidingCenterParallel` """ - doc = rf"""**1. ImplicitDiffusion:** - -{GyrokineticPoissonSolve.__doc__} + doc = rf"""The Poisson solve sits between the two pushers, which both move the guiding centers. With Strang +splitting, bxe(dt/2), poisson, para(dt), poisson, bxe(dt/2), the parallel push and the closing bxe(dt/2) use a freshly +solved potential; the two bxe(dt/2) at a step boundary act as one bxe(dt) with the potential solved before it. With +Lie-Trotter splitting the parallel push uses a fresh potential and the E* x b push the one of the previous step. +Before the scalars and the field data are saved, the Poisson equation is solved for the current marker positions. -**2. PushGuidingCenterBxEstar:** +**1. PushGuidingCenterBxEstar:** {PushGuidingCenterBxEstar.__doc__} +**2. GyrokineticPoissonSolve:** + +{GyrokineticPoissonSolve.__doc__} + **3. PushGuidingCenterParallel:** {PushGuidingCenterParallel.__doc__} From 66bed94ffdd88ee1ddc4d83b4e39bc35a832fc52 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Tue, 20 Oct 2026 08:01:21 +0200 Subject: [PATCH 2/2] DriftKineticElectrostaticAdiabatic: no extra Poisson solve for saving data Remove the update_scalar_quantities override from the previous commit (review of #799): solving the Poisson equation once per saved step is too expensive. The saved potential and en_phi belong to the marker positions before the last push of a time step; this is stated in doc_discretization. Co-Authored-By: Claude Opus 5.5 --- .../drift_kinetic_electrostatic_adiabatic.py | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/src/struphy/models/drift_kinetic_electrostatic_adiabatic.py b/src/struphy/models/drift_kinetic_electrostatic_adiabatic.py index ff70059c2..4f2e11ba3 100644 --- a/src/struphy/models/drift_kinetic_electrostatic_adiabatic.py +++ b/src/struphy/models/drift_kinetic_electrostatic_adiabatic.py @@ -186,18 +186,6 @@ def post_allocate(self): average_matrix_teta = AverageOperator(derham, direction=1) self.average_matrix = average_matrix_z @ average_matrix_teta - def update_scalar_quantities(self): - """Solve the Poisson equation for the current marker positions, then update the scalars. - - A time step ends with a push (bxe with Strang, para with Lie-Trotter), so the potential would lag behind the - markers in ``en_phi``, ``en_tot`` and the saved data (which is extracted after this call). The extra solve, - once per saved step, does not change the dynamics: the next Poisson solve of the time stepping overwrites it. - - :meta private: - """ - self.propagators.gc_poisson(1.0) - super().update_scalar_quantities() - def _compute_en_phi(self): phi = self.em_fields.phi.spline.vector epsilon = self.kinetic_ions.equation_params.epsilon @@ -300,7 +288,8 @@ def doc_discretization(cls): splitting, bxe(dt/2), poisson, para(dt), poisson, bxe(dt/2), the parallel push and the closing bxe(dt/2) use a freshly solved potential; the two bxe(dt/2) at a step boundary act as one bxe(dt) with the potential solved before it. With Lie-Trotter splitting the parallel push uses a fresh potential and the E* x b push the one of the previous step. -Before the scalars and the field data are saved, the Poisson equation is solved for the current marker positions. +Since a time step ends with a push, the saved potential and ``en_phi`` belong to the marker positions before that +last push. **1. PushGuidingCenterBxEstar:**