From 80a06bbb87b1c9c4dcf7e024a9e702e245578f2c Mon Sep 17 00:00:00 2001 From: emile grivet Date: Tue, 2 Jun 2026 16:38:51 +0200 Subject: [PATCH 01/13] first modifications (proper pproc file) --- .../diocotron_instability/params_diocotron.py | 20 +-- .../diocotron_instability/pproc_diocotron.py | 135 ++++++++++-------- 2 files changed, 83 insertions(+), 72 deletions(-) diff --git a/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py b/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py index 1eddd22d4..ddc6d8426 100644 --- a/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py +++ b/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py @@ -69,18 +69,18 @@ ) # List all variables and decide whether to save their data -model.em_fields.phi.save_data = True -model.kinetic_ions.var.save_data = True +model.em_fields.phi.save_data = False +model.kinetic_ions.var.save_data = False # -------------------------- # Instance of the simulation # -------------------------- # Environment options -env = EnvironmentOptions(sim_folder="simdata",profiling_activated=True, profiling_trace=True) +env = EnvironmentOptions(sim_folder="sim_15",profiling_activated=True, profiling_trace=True, restart=True) # Time stepping -time_opts = Time(dt=0.05, Tend=5.2, split_algo="LieTrotter") +time_opts = Time(dt=0.01, Tend=51., split_algo="LieTrotter") # Geometry domain = domains.HollowCylinder(a1=1.0, a2=10.0, Lz=10.0) @@ -116,15 +116,15 @@ # Particle parameters # ------------------- -Np=200000 +Np=500000 loading_params = LoadingParameters(Np = Np, loading="sobol_standard", spatial="disc") -weights_params = WeightsParameters(control_variate=True, reject_weights=True, threshold=0.00001) +weights_params = WeightsParameters(control_variate=True, reject_weights=True, threshold=0.0001) boundary_params = BoundaryParameters() -sorting_params = SortingParameters(boxes_per_dim=(24,24,1), do_sort=True) +sorting_params = SortingParameters(boxes_per_dim=(12,12,1), do_sort=True, sorting_frequency=5) # density binning eta_bin = BinningPlot(slice='e1_e2', n_bins= (128,128), ranges= ((0.0, 1.0), (0.0,1.0))) -saving_params = SavingParameters(binning_plots=(eta_bin, )) +saving_params = SavingParameters(binning_plots=()) model.kinetic_ions.set_markers(loading_params=loading_params, weights_params=weights_params, @@ -139,7 +139,7 @@ # ------------------ model.propagators.gc_poisson.options = model.propagators.gc_poisson.Options() -model.propagators.push_gc_bxe.options = model.propagators.push_gc_bxe.Options(algo="explicit", phi=model.em_fields.phi, evaluate_e_field=True) +model.propagators.push_gc_bxe.options = model.propagators.push_gc_bxe.Options(algo="discrete_gradient_1st_order_newton", phi=model.em_fields.phi, evaluate_e_field=True) # ------------------ # Initial conditions @@ -173,7 +173,7 @@ def n_init(etas,r_minus=r_minus,r_plus=r_plus): # Perturbations for (some) kinetic species # for linear case amps = (1e-6,) -perturbation = perturbations.ModesCos(amps=(0.5,), ms=(ms,), perb_domain=((eta_minus,eta_plus), None, None)) +perturbation = perturbations.ModesCos(amps=(1e-6,), ms=(ms,), perb_domain=((eta_minus,eta_plus), None, None)) init = maxwellians.GyroMaxwellian2D(n=(n_init, perturbation), equil=equil) model.kinetic_ions.var.add_initial_condition(init) diff --git a/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py b/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py index 96c031330..4283bf748 100644 --- a/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py +++ b/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py @@ -1,76 +1,90 @@ -import params_diocotron as params +import importlib.util +import pyvista as pv from struphy import PlottingData, PostProcessor -import os +import os, sys import cunumpy as xp +import scipy.optimize as sc from matplotlib import pyplot as plt import h5py +import logging +from struphy import set_logging_level +set_logging_level(logging.INFO) + # ------------------ # Post process simulation data +# In order to compare different simulations, execute this file as `python pproc_diocotron.py sim_1 sim_2 ...` +# where `sim_1`, `sim_2`, etc. are the names of the simulation folders to be post-processed and plotted together. +# If only one argument, the 2D plots will be shown. If multiple arguments, only the growth rate plot will be shown. # ------------------ def main(): - sim_name = "simdata" - sim_path = os.path.join(os.getcwd(), sim_name) - - pp = PostProcessor(sim=params.sim) - pp.process(physical=True) - - pdata = PlottingData(sim=params.sim) - pdata.load() - - # path to save plots - # save_path = os.path.join(os.getcwd(), "images", "sim") - # os.makedirs(save_path, exist_ok=True) - - # ------------------ - # Check simulation domain - # ------------------ - - params.domain.show() - - # ------------------ - # Determine electrical potentail growth rate - # ------------------ - - # get scalar data (post processing not needed for scalar data) - pa_data = os.path.join(sim_path, "data") - with h5py.File(os.path.join(pa_data, "data_proc0.hdf5"), "r") as f: - time = f["time"]["value"][()] - en_phi = f["scalar"]["en_phi"][()] - - # determine growth rate - exp_func = lambda x,m,b: 10**(m*x+b) - - # time interval to determine growth rate - ti = pdata.t_grid[-1]//4 - if ti == 0.0: - tf = pdata.t_grid[-1] + if len(sys.argv)>1: + sim_names = sys.argv[1:] else: - tf = 2*ti - print(f"{ti = }, {tf = }") - #ti, tf = 2.5, 5.1 - - xi = xp.abs(pdata.t_grid - ti).argmin() + 1 # index of time 100 [a.lu.] (observed end of growth rate) - xf = xp.abs(pdata.t_grid - tf).argmin() + 1 # index of time 200 [a.lu.] (observed end of growth rate) - phi_init=en_phi[1] - en_phi = en_phi - phi_init - fitting = xp.polyfit(time[xi:xf], xp.log10(en_phi[xi:xf]), deg=1) + sim_names = ["sim_5"] + en_phis = [] + times = [] + sls = [] + params_opts = [] + for i, sim_name in enumerate(sim_names): + sim_path = os.path.join(os.getcwd(), sim_name) + + spec = importlib.util.spec_from_file_location("params", os.path.join(sim_path, "parameters.py")) + params = importlib.util.module_from_spec(spec) + spec.loader.exec_module(params) + + if not os.path.isdir(os.path.join(sim_path, "post_processing")): + pp = PostProcessor(sim=params.sim) + pp.process(physical=True) + + pdata = PlottingData(sim=params.sim) + pdata.load() + + # ------------------ + # Determine electrical potentail growth rate + # ------------------ + + # get scalar data (post processing not needed for scalar data) + pa_data = os.path.join(sim_path, "data") + with h5py.File(os.path.join(pa_data, "data_proc0.hdf5"), "r") as f: + times.append(f["time"]["value"][()]) + en_phis.append(xp.power(f["scalar"]["en_phi"][()], 1.0)) + + # time interval to determine growth rate + ti, tf = 0.0, 42.0 + if tf>times[i][-1]: tf = times[i][-1] + if ti>tf: + ti = tf/2 + xi = xp.abs(pdata.t_grid - ti).argmin() # index of time 100 [a.lu.] (observed end of growth rate) + xf = xp.abs(pdata.t_grid - tf).argmin() + 1 # index of time 200 [a.lu.] (observed end of growth rate) + if xi==0: + xi=1 # avoid including t=0 in fit + + sls.append(tuple([slice(xi, xf)])) + + # determine growth rate + fitting_func = lambda x,m,b,c0: xp.exp(m*x+b)+c0 + jac_func = lambda x,m,b,c0: xp.array([x*xp.exp(m*x+b), xp.exp(m*x+b), xp.ones_like(x)]).transpose() + + params_opt, _ = sc.curve_fit(fitting_func, times[i][sls[i]], en_phis[i][sls[i]], p0=(1e-3, -5, en_phis[i][1]), jac=jac_func, maxfev=10000)#3.07e2 + params_opts.append(params_opt) + + logging.info(f"Fitted growth rate for {sim_name}: {params_opt[0]:.4e}") fig, ax = plt.subplots(1, figsize = (18, 12)) - - # plot - ax.plot(time, en_phi, label=r"$\phi$") - ax.plot( - pdata.t_grid, - exp_func(pdata.t_grid, *fitting), - label=f"fitted growth rate {ti=}, {tf=}, growth_rate={fitting[0]:.4e}" - ) + for i in range(len(sim_names)): + ax.scatter(times[i][1:], en_phis[i][1:], label=r"$\phi_{"+sim_names[i][4:]+r"}$", marker='x', s=0.05) + ax.plot( + times[i][sls[i]], + fitting_func(times[i][sls[i]], *params_opts[i]), + label=f"fitted growth rate {ti=}, {tf=}, growth_rate={params_opts[i][0]:.4e}, b={params_opts[i][1]:.4e}, c0={params_opts[i][2]:.4e}" + ) ax.axvline(ti, color="gray", linestyle="--", alpha=0.5) ax.axvline(tf, color="gray", linestyle="--", alpha=0.5) - ax.set_yscale('log') + #ax.set_yscale('log') ax.legend() ax.set_title(f"{params.time_opts.dt=}, {params.time_opts.split_algo=}, {params.grid.num_elements=}, {params.derham_opts.degree=}, {params.loading_params.ppc=}") @@ -79,11 +93,8 @@ def main(): plt.tight_layout() plt.show() - # plt.savefig(os.path.join(save_path, "growth_rate.png")) - # plt.close() - - en_phi = en_phi + phi_init - + if len(sim_names)>1: + exit() # ------------------ # Show evolution of mass density distribution # ------------------ @@ -212,7 +223,7 @@ def extract_images(bin_name, quantity, img_dir): plt.close(fig) # extract_images("e1_e2_density", "f_binned", os.path.join(save_path, "video")) - save_video_pngs = True + save_video_pngs = False if save_video_pngs: if not os.path.exists(sim_path+"/video"): os.mkdir(sim_path+"/video") From 2ffd02797eb1d8882fd94edd2ad88db5ced4a782 Mon Sep 17 00:00:00 2001 From: emile grivet Date: Tue, 2 Jun 2026 16:45:20 +0200 Subject: [PATCH 02/13] proper growth rate obtained --- .../diocotron_instability/params_diocotron.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py b/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py index ddc6d8426..87b750d39 100644 --- a/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py +++ b/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py @@ -77,10 +77,10 @@ # -------------------------- # Environment options -env = EnvironmentOptions(sim_folder="sim_15",profiling_activated=True, profiling_trace=True, restart=True) +env = EnvironmentOptions(sim_folder="sim_1",profiling_activated=True, profiling_trace=True, restart=False) # Time stepping -time_opts = Time(dt=0.01, Tend=51., split_algo="LieTrotter") +time_opts = Time(dt=0.01, Tend=51.0, split_algo="LieTrotter") # Geometry domain = domains.HollowCylinder(a1=1.0, a2=10.0, Lz=10.0) @@ -124,7 +124,7 @@ # density binning eta_bin = BinningPlot(slice='e1_e2', n_bins= (128,128), ranges= ((0.0, 1.0), (0.0,1.0))) -saving_params = SavingParameters(binning_plots=()) +saving_params = SavingParameters(binning_plots=()) # (binning_plots=(eta_bin,)) if you want to save the density binning data model.kinetic_ions.set_markers(loading_params=loading_params, weights_params=weights_params, @@ -172,7 +172,7 @@ def n_init(etas,r_minus=r_minus,r_plus=r_plus): # Perturbations for (some) kinetic species -# for linear case amps = (1e-6,) +# for non linear case amps = (0.5,) perturbation = perturbations.ModesCos(amps=(1e-6,), ms=(ms,), perb_domain=((eta_minus,eta_plus), None, None)) init = maxwellians.GyroMaxwellian2D(n=(n_init, perturbation), equil=equil) model.kinetic_ions.var.add_initial_condition(init) From 37ef6dacab708838164ce31a8938059d42e4e86c Mon Sep 17 00:00:00 2001 From: emile grivet Date: Tue, 2 Jun 2026 16:52:45 +0200 Subject: [PATCH 03/13] formating --- .../ToyGyrokinetic/diocotron_instability/pproc_diocotron.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py b/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py index 4283bf748..68f75a556 100644 --- a/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py +++ b/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py @@ -1,8 +1,8 @@ import importlib.util -import pyvista as pv from struphy import PlottingData, PostProcessor -import os, sys +import os +import sys import cunumpy as xp import scipy.optimize as sc from matplotlib import pyplot as plt From fa6dea219f6e880ff3e1a087df683b8bcc2a04e3 Mon Sep 17 00:00:00 2001 From: emile grivet Date: Wed, 10 Jun 2026 09:41:54 +0200 Subject: [PATCH 04/13] modification of equilibrium --- src/struphy/fields_background/equils.py | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/src/struphy/fields_background/equils.py b/src/struphy/fields_background/equils.py index d047645a7..1b617e065 100644 --- a/src/struphy/fields_background/equils.py +++ b/src/struphy/fields_background/equils.py @@ -932,7 +932,7 @@ class AdhocTorus(AxisymmMHDequilibrium): a : 1. # minor radius R0 : 3. # major radius B0 : 2. # on-axis toroidal magnetic field - q_kind : 0 # which profile (0 : parabolic, 1 : other, see documentation) + q_kind : 0 # which profile (0 : parabolic, 1 : other, 2 : parabolic with linear term, see documentation) q0 : 1.05 # safety factor at r=0 q1 : 1.80 # safety factor at r=a n1 : .5 # 1st shape factor for number density profile @@ -955,6 +955,7 @@ def __init__( q_kind: int = 0, q0: float = 1.71, q1: float = 1.87, + l: float = 0.0, n1: float = 2.0, n2: float = 1.0, na: float = 0.2, @@ -983,7 +984,7 @@ def __init__( self._psi_i = None self._p_i = None - else: + elif self.params["q_kind"] == 1 or self.params["q_kind"] == 2: r_i = xp.linspace(0.0, self.params["a"], self.params["psi_nel"] + 1) def dpsi_dr(r): @@ -1088,7 +1089,7 @@ def psi_r(self, r, der=0): out = self.params["B0"] * (q_bar_0 - r * q_bar_1) / q_bar_0**2 # alternative profile (interpolated) - else: + elif self.params["q_kind"] == 1 or self.params["q_kind"] == 2: out = self._psi_i(r, nu=der) # remove all "dimensions" for point-wise evaluation @@ -1105,6 +1106,7 @@ def q_r(self, r, der=0): q0 = self.params["q0"] q1 = self.params["q1"] + l = self.params["l"] a = self.params["a"] @@ -1114,9 +1116,15 @@ def q_r(self, r, der=0): qout = q0 + (q1 - q0) * (r / a) ** 2 else: qout = 2 * (q1 - q0) * r / a**2 + + elif self.params["q_kind"] == 2: + if der == 0: + qout = q0 + (q1 - q0) * (r / a) ** 2 + l * r / a + else: + qout = 2 * (q1 - q0) * r / a**2 + l / a # alternative profile - else: + elif self.params["q_kind"] == 1: # int/float input if isinstance(r, (int, float)): if r == 0: @@ -1210,7 +1218,7 @@ def p_r(self, r): ) # alternative profile - else: + elif self.params["q_kind"] == 1: pout = self._p_i(r) # remove all "dimensions" for point-wise evaluation @@ -1219,7 +1227,7 @@ def p_r(self, r): pout = pout.item() # ad-hoc profile - else: + elif self.params["p_kind"] == 1: pout = ( self.params["B0"] ** 2 * self.params["beta"] From 637309bbf11a9fb8feddbadb86a27e72fab3df14 Mon Sep 17 00:00:00 2001 From: emile grivet Date: Wed, 10 Jun 2026 09:50:48 +0200 Subject: [PATCH 05/13] add of propagator file --- .../poisson_adiabatic_gyrokinetic.py | 168 ++++++++++++++++++ 1 file changed, 168 insertions(+) create mode 100644 src/struphy/propagators/poisson_adiabatic_gyrokinetic.py diff --git a/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py b/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py new file mode 100644 index 000000000..9fc83cac8 --- /dev/null +++ b/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py @@ -0,0 +1,168 @@ +from dataclasses import dataclass +from typing import Callable, Literal + +from feectools.linalg.stencil import StencilVector + +from struphy.feec.mass import AverageOperator +from struphy.io.options import LiteralOptions +from struphy.linear_algebra.solver import SolverParameters +from struphy.models.variables import FEECVariable +from struphy.pic.accumulation.particles_to_grid import AccumulatorVector +from struphy.pic.base import Particles +from struphy.propagators.implicit_diffusion import ImplicitDiffusion +from struphy.utils.utils import check_option + + +class PoissonAdiabaticGyrokinetic(ImplicitDiffusion): + r""" + Weak discretization of the Poisson equation, with adabatic response of electrons. + + Find :math:`\phi \in H^1` such that + + .. math:: + + \frac{1}{Z\epsilon^2} \int_\Omega \frac{n_0}{T_0} \psi\, \phi\,\textrm d \mathbf x + \int_\Omega \frac{n_0}{|B_0|²} \nabla \psi^\top \, \nabla \phi \,\textrm d \mathbf x = \sum_i \int_\Omega \psi\, \rho_i(\mathbf x)\,\textrm d \mathbf x \qquad \forall \ \psi \in H^1\,, + + where :math:`\epsilon \in \mathbb R` is the gyrokinetic ratio defined in units, and Z the charge number of ions. + Boundary terms from integration by parts are assumed to vanish. + + The equation is discretized as + + .. math:: + + \left( \frac{1}{Z\epsilon^2}\,\mathbb M^0_ad + \mathbb G^\top \mathbb M^1 \mathbb G \right)\, \boldsymbol\phi^{n+1} = \sum_i(\Lambda^0, \rho_i )_{L^2}\,, + + where :math:`\mathbb M^1` is the :math:`H(\textnormal{curl})`-mass matrix + and :math:`\mathbb S` is a stabilization matrix. + + Parameters + ---------- + phi : StencilVector + FE coefficients of the solution as a discrete 0-form. + + stab_eps : float + Stabilization parameter multiplied on stab_mat (default=0.0). + + stab_mat : str + Name of the stabilizing matrix. + + rho : StencilVector or tuple or list + (List of) right-hand side FE coefficients of a 0-form (optional, can be set with a setter later). + Can be either a) StencilVector or b) 2-tuple, or a list of those. + In case b) the first tuple entry must be :class:`~struphy.pic.accumulation.particles_to_grid.AccumulatorVector`, + and the second entry must be :class:`~struphy.pic.base.Particles`. + + x0 : StencilVector + Initial guess for the iterative solver (optional, can be set with a setter later). + + solver : dict + Parameters for the iterative solver (see ``__init__`` for details). + """ + + @dataclass + class Options: + """Configuration options for :class:`PoissonAdiabaticGyrokinetic`. + + Parameters + ---------- + stab_eps : float, default=0.0 + Stabilization weight used for the mass-like term in the Poisson + operator. + Internally mapped to ``sigma_1 = stab_eps`` in the parent + :class:`ImplicitDiffusion` formulation. + + stab_mat : {"M0", "M0ad", "Id"}, default="Id" + Stabilization matrix multiplied by ``stab_eps``. + + - ``"M0"``: standard weighted 0-form mass operator. + - ``"M0ad"``: adiabatic-electron weighted 0-form mass operator. + - ``"Id"``: identity operator. + + rho : FEECVariable or Callable or tuple or list, default=None + Right-hand side source term(s) of the Poisson problem. + Accepted entries are: + + - ``None``: zero source. + - ``FEECVariable`` in ``H1``. + - ``Callable`` to be projected to ``H1`` via ``L2Projector``. + - ``AccumulatorVector``. + - a ``list`` containing any mix of the entries above. + + The tuple form is accepted by typing for compatibility with other + propagator interfaces that pair particle data with accumulators. + + rho_coeffs : float or list, default=None + Multiplicative coefficient(s) applied to ``rho``. + If ``None``, coefficients default to ``1.0`` for all sources. + + x0 : StencilVector, default=None + Initial guess for the iterative linear solver. + + solver : LiteralOptions.OptsSymmSolver, default="pcg" + Name of the symmetric iterative solver passed to + :func:`psydac.linalg.solvers.inverse`. + + precond : LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner" + Name of the preconditioner configuration. + Currently this class inherits the same behavior as + :class:`ImplicitDiffusion`, where ``pc=None`` is used internally. + + solver_params : SolverParameters, default=None + Iterative-solver controls (for example ``tol``, ``maxiter``, + ``verbose``, ``info``, ``recycle``). + If ``None``, defaults to ``SolverParameters()``. + + Notes + ----- + ``Poisson.Options`` reuses :class:`ImplicitDiffusion` internals by + enforcing + ``sigma_2 = 0.0``, ``sigma_3 = 1.0``, ``divide_by_dt = False`` and + ``diffusion_mat = "M1"`` in ``__post_init__``. + """ + + # specific literals + OptsStabMat = Literal["M0", "M0ad", "Id", "M0ad_withT"] + OptsDiffusionMat = Literal["M1", "M1perp", "M1gyro"] + # propagator options + epsilon: float = 0.0 + stab_mat: OptsStabMat = "M0ad_withT" + diffusion_mat = "M1gyro" + rho: FEECVariable | Callable | tuple[AccumulatorVector, Particles] | list = None + rho_coeffs: float | list = None + x0: StencilVector = None + solver: LiteralOptions.OptsSymmSolver = "pcg" + precond: LiteralOptions.OptsMassPrecond = "MassMatrixPreconditioner" + solver_params: SolverParameters = None + + def __post_init__(self): + # checks + check_option(self.stab_mat, self.OptsStabMat) + check_option(self.solver, LiteralOptions.OptsSymmSolver) + check_option(self.precond, LiteralOptions.OptsMassPrecond) + + # defaults + if self.solver_params is None: + self.solver_params = SolverParameters() + + # Poisson solve (-> set some params of parent class) + self.sigma_1 = 1.0 + self.sigma_2 = 0.0 + self.sigma_3 = 1.0 + self.divide_by_dt = False + + def allocate(self, verbose = False): + super().allocate(verbose) + average_mat = AverageOperator(self.derham, "H1", 2) + temp = self._stab_mat.copy() @ average_mat + self._stab_mat -= temp + + @property + def options(self) -> Options: + if not hasattr(self, "_options"): + self._options = self.Options() + return self._options + + @options.setter + def options(self, new): + assert isinstance(new, self.Options) + self._options = new From abd8476a2f069fe1715fa4f3475260976755400b Mon Sep 17 00:00:00 2001 From: emile grivet Date: Wed, 10 Jun 2026 09:57:31 +0200 Subject: [PATCH 06/13] docstring equil --- src/struphy/fields_background/equils.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/struphy/fields_background/equils.py b/src/struphy/fields_background/equils.py index 1b617e065..b51363196 100644 --- a/src/struphy/fields_background/equils.py +++ b/src/struphy/fields_background/equils.py @@ -871,6 +871,9 @@ class AdhocTorus(AxisymmMHDequilibrium): &q_0 + ( q_1 - q_0 )\frac{r^2}{a^2} \quad &&\textnormal{if} \quad q_\textnormal{kind}=0\,, &\frac{q_0}{1-\left(1-\frac{r^2}{a^2}\right)^{\frac{q_1}{q_0}}}\frac{r^2}{a^2} \quad &&\textnormal{if} \quad q_\textnormal{kind}=1\,. + + + &q_0 + l\frac{r}{a} ( q_1 - q_0 )\frac{r^2}{a^2} \quad &&\textnormal{if} \quad q_\textnormal{kind}=2\,, \end{aligned}\right. The pressure profile @@ -903,6 +906,8 @@ class AdhocTorus(AxisymmMHDequilibrium): Safety factor at r=0 (default: 1.71). q1 : float Safety factor at r=a (default: 1.87). + l : float + Linear term factor for q profile if q_kind=2 (default: 0.). n1 : float 1st shape factor for ion number density profile (default: 0.). n2 : float @@ -935,6 +940,7 @@ class AdhocTorus(AxisymmMHDequilibrium): q_kind : 0 # which profile (0 : parabolic, 1 : other, 2 : parabolic with linear term, see documentation) q0 : 1.05 # safety factor at r=0 q1 : 1.80 # safety factor at r=a + l : 0. # linear term factor for q profile if q_kind=2 n1 : .5 # 1st shape factor for number density profile n2 : 1. # 2nd shape factor for number density profile na : .2 # number density at r=a From 945d6855ca2316ffa5b5098dcc42735e5977fad0 Mon Sep 17 00:00:00 2001 From: emile grivet Date: Wed, 10 Jun 2026 09:58:38 +0200 Subject: [PATCH 07/13] formating --- src/struphy/fields_background/equils.py | 2 +- src/struphy/propagators/poisson_adiabatic_gyrokinetic.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/struphy/fields_background/equils.py b/src/struphy/fields_background/equils.py index b51363196..aed6175f2 100644 --- a/src/struphy/fields_background/equils.py +++ b/src/struphy/fields_background/equils.py @@ -1122,7 +1122,7 @@ def q_r(self, r, der=0): qout = q0 + (q1 - q0) * (r / a) ** 2 else: qout = 2 * (q1 - q0) * r / a**2 - + elif self.params["q_kind"] == 2: if der == 0: qout = q0 + (q1 - q0) * (r / a) ** 2 + l * r / a diff --git a/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py b/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py index 9fc83cac8..4bee87b50 100644 --- a/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py +++ b/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py @@ -150,7 +150,7 @@ def __post_init__(self): self.sigma_3 = 1.0 self.divide_by_dt = False - def allocate(self, verbose = False): + def allocate(self, verbose=False): super().allocate(verbose) average_mat = AverageOperator(self.derham, "H1", 2) temp = self._stab_mat.copy() @ average_mat From 94cdafa217173ed58081a7f2e37f462d37eadfda Mon Sep 17 00:00:00 2001 From: emile grivet Date: Wed, 10 Jun 2026 10:12:55 +0200 Subject: [PATCH 08/13] adddition of new propagator in doc --- doc/sections/subsections/propagators-avail.rst | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/doc/sections/subsections/propagators-avail.rst b/doc/sections/subsections/propagators-avail.rst index b2227b52b..bbdd19c28 100644 --- a/doc/sections/subsections/propagators-avail.rst +++ b/doc/sections/subsections/propagators-avail.rst @@ -11,6 +11,11 @@ Field solvers :exclude-members: options, allocate :show-inheritance: +.. automodule:: struphy.propagators.poisson_adiabatic_gyrokinetic + :members: + :exclude-members: options, allocate + :show-inheritance: + .. automodule:: struphy.propagators.curl_curl_solve :members: :exclude-members: options, allocate From 81df61b162fb76af880c66d63ec30feb1e751c2c Mon Sep 17 00:00:00 2001 From: emile grivet Date: Wed, 10 Jun 2026 10:23:04 +0200 Subject: [PATCH 09/13] diocotron example postprocessing rectification --- .../diocotron_instability/pproc_diocotron.py | 2 +- .../propagators/poisson_adiabatic_gyrokinetic.py | 14 +++++++++++--- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py b/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py index 68f75a556..49f8cd8fb 100644 --- a/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py +++ b/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py @@ -23,7 +23,7 @@ def main(): if len(sys.argv)>1: sim_names = sys.argv[1:] else: - sim_names = ["sim_5"] + sim_names = ["sim_1"] en_phis = [] times = [] sls = [] diff --git a/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py b/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py index 4bee87b50..f3648f9f4 100644 --- a/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py +++ b/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py @@ -123,9 +123,11 @@ class Options: # specific literals OptsStabMat = Literal["M0", "M0ad", "Id", "M0ad_withT"] OptsDiffusionMat = Literal["M1", "M1perp", "M1gyro"] + OptsGeometry = Literal["cylindrical", "toroidal"] # propagator options epsilon: float = 0.0 stab_mat: OptsStabMat = "M0ad_withT" + which_geometry: OptsGeometry = "cylindrical" diffusion_mat = "M1gyro" rho: FEECVariable | Callable | tuple[AccumulatorVector, Particles] | list = None rho_coeffs: float | list = None @@ -152,9 +154,15 @@ def __post_init__(self): def allocate(self, verbose=False): super().allocate(verbose) - average_mat = AverageOperator(self.derham, "H1", 2) - temp = self._stab_mat.copy() @ average_mat - self._stab_mat -= temp + if self.option.which_geometry == "cylindrical": + average_mat = AverageOperator(self.derham, "H1", 2) + temp = self._stab_mat.copy() @ average_mat + self._stab_mat -= temp + elif self.option.which_geometry == "toroidal": + average_mat1 = AverageOperator(self.derham, "H1", 1) + average_mat2 = AverageOperator(self.derham, "H1", 2) + temp = self._stab_mat.copy() @ average_mat1 @ average_mat2 + self._stab_mat -= temp @property def options(self) -> Options: From 59e4b1f84dcdfb687d00f207edeb3add12ddba0f Mon Sep 17 00:00:00 2001 From: emile grivet Date: Wed, 10 Jun 2026 14:21:26 +0200 Subject: [PATCH 10/13] rectification pproc_diocotron --- .../diocotron_instability/pproc_diocotron.py | 38 ++++++++++++------- 1 file changed, 24 insertions(+), 14 deletions(-) diff --git a/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py b/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py index 49f8cd8fb..1865b2aaa 100644 --- a/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py +++ b/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py @@ -10,7 +10,7 @@ import logging from struphy import set_logging_level -set_logging_level(logging.INFO) +set_logging_level(logging.DEBUG) # ------------------ @@ -20,14 +20,18 @@ # If only one argument, the 2D plots will be shown. If multiple arguments, only the growth rate plot will be shown. # ------------------ def main(): - if len(sys.argv)>1: - sim_names = sys.argv[1:] + if __name__ == "main": + if len(sys.argv)>1: + sim_names = sys.argv[1:] + else: + sim_names = ["sim_1"] else: sim_names = ["sim_1"] en_phis = [] times = [] sls = [] params_opts = [] + fitting = [] for i, sim_name in enumerate(sim_names): sim_path = os.path.join(os.getcwd(), sim_name) @@ -42,6 +46,7 @@ def main(): pdata = PlottingData(sim=params.sim) pdata.load() + logging.warning(f"arriving at post processing with these sim_names : {sim_names}, knowing that {__name__=}") # ------------------ # Determine electrical potentail growth rate # ------------------ @@ -64,23 +69,28 @@ def main(): sls.append(tuple([slice(xi, xf)])) - # determine growth rate - fitting_func = lambda x,m,b,c0: xp.exp(m*x+b)+c0 - jac_func = lambda x,m,b,c0: xp.array([x*xp.exp(m*x+b), xp.exp(m*x+b), xp.ones_like(x)]).transpose() + if len(times) > 3: + fitting.append(True) + # determine growth rate + fitting_func = lambda x,m,b,c0: xp.exp(m*x+b)+c0 + jac_func = lambda x,m,b,c0: xp.array([x*xp.exp(m*x+b), xp.exp(m*x+b), xp.ones_like(x)]).transpose() - params_opt, _ = sc.curve_fit(fitting_func, times[i][sls[i]], en_phis[i][sls[i]], p0=(1e-3, -5, en_phis[i][1]), jac=jac_func, maxfev=10000)#3.07e2 - params_opts.append(params_opt) + params_opt, _ = sc.curve_fit(fitting_func, times[i][sls[i]], en_phis[i][sls[i]], p0=(1e-3, -5, en_phis[i][1]), jac=jac_func, maxfev=10000)#3.07e2 + params_opts.append(params_opt) - logging.info(f"Fitted growth rate for {sim_name}: {params_opt[0]:.4e}") + logging.info(f"Fitted growth rate for {sim_name}: {params_opt[0]:.4e}") + else: + fitting.append(False) fig, ax = plt.subplots(1, figsize = (18, 12)) for i in range(len(sim_names)): ax.scatter(times[i][1:], en_phis[i][1:], label=r"$\phi_{"+sim_names[i][4:]+r"}$", marker='x', s=0.05) - ax.plot( - times[i][sls[i]], - fitting_func(times[i][sls[i]], *params_opts[i]), - label=f"fitted growth rate {ti=}, {tf=}, growth_rate={params_opts[i][0]:.4e}, b={params_opts[i][1]:.4e}, c0={params_opts[i][2]:.4e}" - ) + if fitting[i]: + ax.plot( + times[i][sls[i]], + fitting_func(times[i][sls[i]], *params_opts[i]), + label=f"fitted growth rate {ti=}, {tf=}, growth_rate={params_opts[i][0]:.4e}, b={params_opts[i][1]:.4e}, c0={params_opts[i][2]:.4e}" + ) ax.axvline(ti, color="gray", linestyle="--", alpha=0.5) ax.axvline(tf, color="gray", linestyle="--", alpha=0.5) From 05fdf0a3c267ba1f1cf398513cdaf60d2e1399cc Mon Sep 17 00:00:00 2001 From: emile grivet Date: Wed, 10 Jun 2026 14:24:15 +0200 Subject: [PATCH 11/13] formating --- .../ToyGyrokinetic/diocotron_instability/pproc_diocotron.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py b/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py index 1865b2aaa..002411190 100644 --- a/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py +++ b/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py @@ -10,7 +10,7 @@ import logging from struphy import set_logging_level -set_logging_level(logging.DEBUG) +set_logging_level(logging.INFO) # ------------------ From 853d49ed0abed1dcf861d88ebdc8bed44eeef2c9 Mon Sep 17 00:00:00 2001 From: emile grivet Date: Wed, 10 Jun 2026 14:38:20 +0200 Subject: [PATCH 12/13] rectification pparam_diocotron --- .../ToyGyrokinetic/diocotron_instability/params_diocotron.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py b/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py index 87b750d39..197859ce5 100644 --- a/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py +++ b/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py @@ -69,7 +69,7 @@ ) # List all variables and decide whether to save their data -model.em_fields.phi.save_data = False +model.em_fields.phi.save_data = True model.kinetic_ions.var.save_data = False # -------------------------- @@ -124,7 +124,7 @@ # density binning eta_bin = BinningPlot(slice='e1_e2', n_bins= (128,128), ranges= ((0.0, 1.0), (0.0,1.0))) -saving_params = SavingParameters(binning_plots=()) # (binning_plots=(eta_bin,)) if you want to save the density binning data +saving_params = SavingParameters(binning_plots=(eta_bin,)) # (binning_plots=(eta_bin,)) if you want to save the density binning data model.kinetic_ions.set_markers(loading_params=loading_params, weights_params=weights_params, From 9405b7945cdbb28e22bd4202ff1e21a81c97a1b2 Mon Sep 17 00:00:00 2001 From: emile grivet Date: Thu, 11 Jun 2026 10:34:26 +0200 Subject: [PATCH 13/13] model and poisson prop modified --- .../drift_kinetic_electrostatic_adiabatic.py | 37 ++++++++++++++++--- .../poisson_adiabatic_gyrokinetic.py | 23 ++++++------ 2 files changed, 42 insertions(+), 18 deletions(-) diff --git a/src/struphy/models/drift_kinetic_electrostatic_adiabatic.py b/src/struphy/models/drift_kinetic_electrostatic_adiabatic.py index 2e50838c9..eb82c25f7 100644 --- a/src/struphy/models/drift_kinetic_electrostatic_adiabatic.py +++ b/src/struphy/models/drift_kinetic_electrostatic_adiabatic.py @@ -4,7 +4,7 @@ from feectools.ddm.mpi import mpi as MPI from struphy import BaseUnits -from struphy.feec.mass import L2Projector +from struphy.feec.mass import AverageOperator, L2Projector from struphy.io.options import LiteralOptions from struphy.kinetic_background.base import KineticBackground from struphy.models.base import StruphyModel @@ -17,7 +17,7 @@ from struphy.pic.accumulation import accum_kernels_gc from struphy.pic.accumulation.particles_to_grid import AccumulatorVector from struphy.propagators.base import Propagator -from struphy.propagators.implicit_diffusion import ImplicitDiffusion +from struphy.propagators.poisson_adiabatic_gyrokinetic import PoissonAdiabaticGyrokinetic from struphy.propagators.push_guiding_center_bx_estar import PushGuidingCenterBxEstar from struphy.propagators.push_guiding_center_parallel import PushGuidingCenterParallel from struphy.utils.pyccel import Pyccelkernel @@ -64,19 +64,21 @@ def __init__( charge_number: int = 1, mass_number: float = 1.0, epsilon: float = None, + alpha: float = None, ): self.var = PICVariable(space="Particles5D") self.init_variables( charge_number=charge_number, mass_number=mass_number, epsilon=epsilon, + alpha=alpha, ) ## propagators class Propagators: def __init__(self): - self.gc_poisson = ImplicitDiffusion() + self.gc_poisson = PoissonAdiabaticGyrokinetic() self.push_gc_bxe = PushGuidingCenterBxEstar() self.push_gc_para = PushGuidingCenterParallel() @@ -88,6 +90,7 @@ def __init__( charge_number: int = 1, mass_number: float = 1.0, epsilon: float = None, + alpha: float = None, ): # 0. store input parameters @@ -99,6 +102,7 @@ def __init__( charge_number, mass_number, epsilon, + alpha, ) # 2. derive units (must be done after instantiating species to access charge and mass numbers) @@ -114,10 +118,12 @@ def __init__( # 5. define scalars to be tracked during simulation field_energy = FunctionScalarFEEC(self._compute_en_phi) + phi_integral_ITGtest = FunctionScalarFEEC(self._compute_phi_integral_ITGtest) particle_kinetic = KineticEnergyPIC(self.kinetic_ions.var) particle_magnetic = FunctionScalarPIC(self._compute_en_particle_magnetic, self.kinetic_ions.var) particle_energy = particle_kinetic + particle_magnetic self.scalars = Scalars( + phi_integral=phi_integral_ITGtest, en_phi=field_energy, en_particles=particle_energy, en_tot=field_energy + particle_energy, @@ -170,10 +176,22 @@ def allocate_helpers(self): self.propagators.gc_poisson.options.sigma_1 = 1.0 / epsilon**2 / Z self.propagators.gc_poisson.options.sigma_2 = 0.0 self.propagators.gc_poisson.options.sigma_3 = 1.0 / epsilon - self.propagators.gc_poisson.options.stab_mat = "M0ad" - self.propagators.gc_poisson.options.diffusion_mat = "M1perp" + self.propagators.gc_poisson.options.stab_mat = "M0ad_withT" + self.propagators.gc_poisson.options.diffusion_mat = "M1gyro" self.propagators.gc_poisson.options.rho = rho + self.propagators.gc_poisson.options.divide_by_dt = False self.propagators.gc_poisson.allocate() + self.propagators.gc_poisson(1.0) + + # allocate temporary tabs for scalars computation + derham = self.em_fields.phi.spline.derham + self.phi_squared_averaged_feec = FEECVariable("H1") + self.phi_squared_averaged_feec.allocate(derham, self.em_fields.phi.spline.domain) + self.phi_squared = self.em_fields.phi.spline.vector.space.zeros() + self.phi_squared_average = self.em_fields.phi.spline.vector.space.zeros() + average_matrix_z = AverageOperator(derham, direction=2) + average_matrix_teta = AverageOperator(derham, direction=1) + self.average_matrix = average_matrix_z @ average_matrix_teta def _compute_en_phi(self): phi = self.em_fields.phi.spline.vector @@ -184,6 +202,13 @@ def _compute_en_phi(self): en_phi = 0.5 / epsilon**2 * Propagator.mass_ops.M0ad.dot_inner(phi, phi) return en_phi + en_phi1 + def _compute_phi_integral_ITGtest(self): + phi = self.em_fields.phi.spline.vector + self.phi_squared._data = phi._data**2 + self.average_matrix.dot(self.phi_squared, out=self.phi_squared_average) + self.phi_squared_averaged_feec.spline.vector = self.phi_squared_average + return self.phi_squared_averaged_feec.spline(0.5, 0.5, 0.5)[0, 0, 0] * 2 * xp.pi * 1506.759067 + def _compute_en_particle_magnetic(self): particles = self.kinetic_ions.var.particles particles.save_magnetic_background_energy() @@ -273,7 +298,7 @@ def doc_discretization(cls): """ doc = rf"""**1. ImplicitDiffusion:** -{ImplicitDiffusion.__doc__} +{PoissonAdiabaticGyrokinetic.__doc__} **2. PushGuidingCenterBxEstar:** diff --git a/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py b/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py index f3648f9f4..4bf1fcb4d 100644 --- a/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py +++ b/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py @@ -65,12 +65,6 @@ class Options: Parameters ---------- - stab_eps : float, default=0.0 - Stabilization weight used for the mass-like term in the Poisson - operator. - Internally mapped to ``sigma_1 = stab_eps`` in the parent - :class:`ImplicitDiffusion` formulation. - stab_mat : {"M0", "M0ad", "Id"}, default="Id" Stabilization matrix multiplied by ``stab_eps``. @@ -78,6 +72,12 @@ class Options: - ``"M0ad"``: adiabatic-electron weighted 0-form mass operator. - ``"Id"``: identity operator. + which_geometry: {"cylindrical", "toroidal"}, default="cylindrical" + Geometry of the problem, determines the meaning of `<\phi>` in the stabilization term. + + diffusion_mat : {"M1", "M1perp", "M1gyro"}, defaults="M1gyro" + Diffusion matrix. + rho : FEECVariable or Callable or tuple or list, default=None Right-hand side source term(s) of the Poisson problem. Accepted entries are: @@ -125,10 +125,9 @@ class Options: OptsDiffusionMat = Literal["M1", "M1perp", "M1gyro"] OptsGeometry = Literal["cylindrical", "toroidal"] # propagator options - epsilon: float = 0.0 stab_mat: OptsStabMat = "M0ad_withT" which_geometry: OptsGeometry = "cylindrical" - diffusion_mat = "M1gyro" + diffusion_mat: OptsDiffusionMat = "M1gyro" rho: FEECVariable | Callable | tuple[AccumulatorVector, Particles] | list = None rho_coeffs: float | list = None x0: StencilVector = None @@ -152,13 +151,13 @@ def __post_init__(self): self.sigma_3 = 1.0 self.divide_by_dt = False - def allocate(self, verbose=False): - super().allocate(verbose) - if self.option.which_geometry == "cylindrical": + def allocate(self): + super().allocate() + if self.options.which_geometry == "cylindrical": average_mat = AverageOperator(self.derham, "H1", 2) temp = self._stab_mat.copy() @ average_mat self._stab_mat -= temp - elif self.option.which_geometry == "toroidal": + elif self.options.which_geometry == "toroidal": average_mat1 = AverageOperator(self.derham, "H1", 1) average_mat2 = AverageOperator(self.derham, "H1", 2) temp = self._stab_mat.copy() @ average_mat1 @ average_mat2