[docs]
class CurrentCoupling5DDensity(Propagator):
r""":ref:`FEEC <gempic>` discretization of the following equations:
find :math:`\mathbf U \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}` and :math:`\mathbf B \in H(\textnormal{div})` such that
.. math::
\int \rho_0 \frac{\partial \tilde{\mathbf U}}{\partial t} \cdot \mathbf V \, \textnormal{d} \mathbf{x} = \frac{A_\textnormal{h}}{A_b} \frac{1}{\epsilon} \iiint f^\text{vol} \left(1 - \frac{B_\parallel}{B^*_\parallel}\right) \tilde{\mathbf U} \times \mathbf B_f \cdot \mathbf V \, \textnormal{d} \mathbf{x} \textnormal{d} v_\parallel \textnormal{d} \mu \quad \forall \,\mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\,,
FE coefficients update:
.. math::
\mathbf u^{n+1} - \mathbf u^n = -\frac{A_\textnormal{h}}{A_b} \frac{1}{\epsilon} \mathbb{L}²{\mathbb{B}}^\times_f \mathbb{N}(1/g) \mathbb{W} \mathbb{N}\left(1- \frac{\hat B^0_\parallel}{\hat B^{*0} _\parallel}\right) (\mathbb{L}²)^\top \frac{\Delta t}{2} \cdot (\mathbf u^{n+1} + \mathbf u^n) \,.
For the detail explanation of the notations, see `2022_DriftKineticCurrentCoupling <https://gitlab.mpcdf.mpg.de/struphy/struphy-projects/-/blob/main/running-projects/2022_DriftKineticCurrentCoupling.md?ref_type=heads>`_.
"""
[docs]
class Variables:
"""Container for variables advanced by :class:`CurrentCoupling5DDensity`.
Attributes
----------
u : FEECVariable
Fluid-like FEEC variable in one of ``"Hcurl"``, ``"Hdiv"``, or
``"H1vec"``.
"""
def __init__(self):
self._u: FEECVariable = None
@property
def u(self) -> FEECVariable:
return self._u
@u.setter
def u(self, new):
assert isinstance(new, FEECVariable)
assert new.space in ("Hcurl", "Hdiv", "H1vec")
self._u = new
def __init__(
self,
energetic_ions: PICVariable = None,
b_tilde: FEECVariable = None,
):
"""
Parameters
----------
energetic_ions : PICVariable, default=None
Energetic-ion particle distribution (``"Particles5D"`` space) providing
the density source term.
b_tilde : FEECVariable, default=None
Magnetic perturbation 1-form (``"Hcurl"`` space) entering the coupling term.
If ``None``, only the equilibrium field is used.
"""
self.variables = self.Variables()
self.energetic_ions = energetic_ions
self.b_tilde = b_tilde
[docs]
@dataclass(repr=False)
class Options(OptionsBase):
r"""Configuration options for :class:`CurrentCoupling5DDensity`.
Parameters
----------
ep_scale : float, default=None
Factor of the energetic-particle contributions (accumulation): the mass ratio
:math:`A_\textnormal{h}/A_\textnormal{b}` of the energetic ions and the MHD bulk if ``None``.
u_space : LiteralOptions.OptsVecSpace, default="Hdiv"
FEEC space used for the unknown ``u`` variable.
solver : LiteralOptions.OptsGenSolver, default="pbicgstab"
General iterative solver used for the linear system.
precond : LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner"
Preconditioner applied to the FEEC mass matrix.
solver_params : SolverParameters, default=None
Iterative-solver controls. If ``None``, defaults to
``SolverParameters()``.
filter_params : FilterParameters, default=None
Filter of the accumulated (particle-to-grid) vectors; defaults to ``FilterParameters()`` (no filter).
The three-point (binomial) filter ``FilterParameters(use_filter="three_point")`` reduces the noise,
but the propagator then no longer conserves the energy exactly (the filter is applied to the
accumulated vectors only, not to the fields seen by the particles).
matrix_free : bool, default=True
Apply the particle matrix :math:`\sum_p \vec \Lambda_p F_p \vec \Lambda_p^\top` (antisymmetric 3x3
:math:`F_p`) in each iteration of the solver from data stored per marker
(:class:`~struphy.pic.accumulation.particles_to_grid.MarkerOperator`) instead of assembling it: one spline
evaluation and one deposition per marker and iteration instead of 3 blocks of size
:math:`n_\textnormal{basis}^2` per marker. Same result up to the solver tolerance.
"""
# propagator options
ep_scale: float = None
u_space: LiteralOptions.OptsVecSpace = "Hdiv"
solver: LiteralOptions.OptsGenSolver = "pbicgstab"
precond: LiteralOptions.OptsMassPrecond = "MassMatrixPreconditioner"
solver_params: SolverParameters = None
filter_params: FilterParameters = None
matrix_free: bool = True
def __post_init__(self):
# checks
assert isinstance(self.matrix_free, bool)
check_option(self.u_space, LiteralOptions.OptsVecSpace)
check_option(self.solver, LiteralOptions.OptsGenSolver)
check_option(self.precond, LiteralOptions.OptsMassPrecond)
assert self.ep_scale is None or isinstance(self.ep_scale, float)
# defaults
if self.solver_params is None:
self.solver_params = SolverParameters()
if self.filter_params is None:
self.filter_params = FilterParameters()
@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
logger.info(f"\nNew options for propagator '{self.__class__.__name__}':\n{self._options}")
@profile
def allocate(self):
# factor of the energetic-particle contributions (default: mass ratio A_h/A_b)
self._ep_scale = self.options.ep_scale
if self._ep_scale is None:
self._ep_scale = self.energetic_ions.species.mass_number / self.variables.u.species.mass_number
logger.info(f"{self.__class__.__name__}: ep_scale = A_h/A_b = {self._ep_scale} (mass ratio, default).")
if self.options.u_space == "H1vec":
self._u_form_int = 0
else:
self._u_form_int = int(self.derham.space_to_form[self.options.u_space])
# call operatros
id_M = "M" + self.derham.space_to_form[self.options.u_space] + "n"
self._A = getattr(self.mass_ops, id_M)
# magnetic equilibrium field
unit_b1 = self.projected_equil.unit_b1
curl_unit_b2 = self.projected_equil.curl_unit_b2
self._b2 = self.projected_equil.b2
# scaling factor
epsilon = self.energetic_ions.species.equation_params.epsilon
# temporary vectors to avoid memory allocation
self._b_full = self._b2.space.zeros()
self._rhs_v = self.variables.u.spline.vector.space.zeros()
self._u_new = self.variables.u.spline.vector.space.zeros()
# define Accumulator and arguments
particles = self.energetic_ions.particles
self._args_accum_kernel = (
epsilon,
self._ep_scale,
self._b_full[0]._data,
self._b_full[1]._data,
self._b_full[2]._data,
unit_b1[0]._data,
unit_b1[1]._data,
unit_b1[2]._data,
curl_unit_b2[0]._data,
curl_unit_b2[1]._data,
curl_unit_b2[2]._data,
self._u_form_int,
)
if self.options.matrix_free:
# the 3x3 filling matrices per marker (in the first free columns of the marker array)
column = particles.args_markers.first_free_idx
assert particles.n_cols_aux >= 4, "matrix_free needs 4 free (auxiliary) marker columns."
self._args_accum_kernel += (column,)
# the kernel that stores the filling matrices (called in __call__, kept for kernels())
self._store_kernel = cc_lin_mhd_5d_D_mf
self._A_particles = MarkerOperator(
particles,
self.options.u_space,
marker_operator_asym_dot,
self.mass_ops,
self.domain.args_domain,
(self._u_form_int, column),
symmetry="asym",
)
else:
self._ACC = Accumulator(
particles,
self.options.u_space,
cc_lin_mhd_5d_D,
self.mass_ops,
self.domain.args_domain,
add_vector=False,
symmetry="asym",
filter_params=self.options.filter_params,
)
self._A_particles = self._ACC.operators[0]
# Preconditioner
if self.options.precond is None:
pc = None
else:
pc_class = getattr(preconditioner, self.options.precond)
pc = pc_class(getattr(self.mass_ops, id_M))
# linear solver
self._A_inv = inverse(
self._A,
self.options.solver,
pc=pc,
tol=self.options.solver_params.tol,
maxiter=self.options.solver_params.maxiter,
verbose=self.options.solver_params.verbose,
)
def __call__(self, dt):
# current FE coeffs
un = self.variables.u.spline.vector
# sum up total magnetic field b_full1 = b_eq + b_tilde (in-place)
b_full = self._b2.copy(out=self._b_full)
if self.b_tilde is not None:
b_full += self.b_tilde.spline.vector
b_full.update_ghost_regions()
if self.options.matrix_free:
particles = self.energetic_ions.particles
with ProfileManager.profile_region("kernel: " + self._store_kernel.name):
self._store_kernel(
particles.args_markers,
self.derham.args_derham,
self.domain.args_domain,
*self._args_accum_kernel,
n_threads=particles.args_markers.n_markers,
)
else:
self._ACC(
*self._args_accum_kernel,
)
# define system (M - dt/2 * A)*u^(n + 1) = (M + dt/2 * A)*u^n
lhs = self._A - dt / 2 * self._A_particles
rhs = self._A + dt / 2 * self._A_particles
# solve linear system for updated u coefficients (in-place)
rhs = rhs.dot(un, out=self._rhs_v)
self._A_inv.linop = lhs
with ProfileManager.profile_region(self._solve_region, functions=[self._A_inv.solve]):
_u = self._A_inv.solve(rhs, out=self._u_new)
info = self._A_inv._info
diffs = self.update_feec_variables(u=_u)
if self.options.solver_params.info and MPI.COMM_WORLD.Get_rank() == 0:
logger.info(f"Status for CurrentCoupling5DDensity: {info['success']}")
logger.info(f"Iterations for CurrentCoupling5DDensity: {info['niter']}")
logger.info(f"Maxdiff up for CurrentCoupling5DDensity: {diffs['u']}")
logger.info("")
