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 diff --git a/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py b/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py index 1eddd22d4..197859ce5 100644 --- a/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py +++ b/examples/ToyGyrokinetic/diocotron_instability/params_diocotron.py @@ -70,17 +70,17 @@ # 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.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_1",profiling_activated=True, profiling_trace=True, restart=False) # Time stepping -time_opts = Time(dt=0.05, Tend=5.2, 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) @@ -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=(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, @@ -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 @@ -172,8 +172,8 @@ 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)) +# 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) diff --git a/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py b/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py index 96c031330..002411190 100644 --- a/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py +++ b/examples/ToyGyrokinetic/diocotron_instability/pproc_diocotron.py @@ -1,76 +1,100 @@ -import params_diocotron as params +import importlib.util from struphy import PlottingData, PostProcessor import os +import 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 __name__ == "main": + if len(sys.argv)>1: + sim_names = sys.argv[1:] + else: + sim_names = ["sim_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_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) + + 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() + + logging.warning(f"arriving at post processing with these sim_names : {sim_names}, knowing that {__name__=}") + # ------------------ + # 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)])) + + 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) + + 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)) - - # 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) + 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) - 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 +103,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 +233,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") diff --git a/src/struphy/fields_background/equils.py b/src/struphy/fields_background/equils.py index d047645a7..aed6175f2 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 @@ -932,9 +937,10 @@ 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 + 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 @@ -955,6 +961,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 +990,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 +1095,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 +1112,7 @@ def q_r(self, r, der=0): q0 = self.params["q0"] q1 = self.params["q1"] + l = self.params["l"] a = self.params["a"] @@ -1115,8 +1123,14 @@ def q_r(self, r, der=0): 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 +1224,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 +1233,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"] 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 new file mode 100644 index 000000000..4bf1fcb4d --- /dev/null +++ b/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py @@ -0,0 +1,175 @@ +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_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. + + 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: + + - ``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"] + OptsGeometry = Literal["cylindrical", "toroidal"] + # propagator options + stab_mat: OptsStabMat = "M0ad_withT" + which_geometry: OptsGeometry = "cylindrical" + diffusion_mat: OptsDiffusionMat = "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): + 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.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 + 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