Particle-field coupling propagators

Particle and FEEC variables are updated.

class struphy.propagators.propagators_coupling.VlasovAmpere

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf E \in H(\textnormal{curl})\) and \(f\) such that

\[\begin{split}-& \int_\Omega \frac{\partial \mathbf E}{\partial t} \cdot \mathbf F\,\textrm d \mathbf x = \frac{\alpha^2}{\varepsilon} \int_\Omega \int_{\mathbb{R}^3} f \mathbf{v} \cdot \mathbf F \, \text{d}^3 \mathbf{v} \,\textrm d \mathbf x \qquad \forall \, \mathbf F \in H(\textnormal{curl}) \,, \\[2mm] &\frac{\partial f}{\partial t} + \frac{1}{\varepsilon}\, \mathbf{E} \cdot \frac{\partial f}{\partial \mathbf{v}} = 0 \,.\end{split}\]

Time discretization: Crank-Nicolson (implicit mid-point). System size reduction via SchurSolver, such that

\[\begin{split}\begin{bmatrix} \mathbb{M}_1 \left( \mathbf{e}^{n+1} - \mathbf{e}^n \right) \\ \mathbf{V}^{n+1} - \mathbf{V}^n \end{bmatrix} = \frac{\Delta t}{2} \begin{bmatrix} 0 & - \frac{\alpha^2}{\varepsilon} \mathbb L^1 \bar{DF^{-1}} \bar{\mathbf w} \\ \frac{1}{\varepsilon} \bar{DF^{-\top}} \left(\mathbb L^1\right)^\top & 0 \end{bmatrix} \begin{bmatrix} \mathbf{e}^{n+1} + \mathbf{e}^n \\ \mathbf{V}^{n+1} + \mathbf{V}^n \end{bmatrix}\end{split}\]

based on the SchurSolver with

\[A = \mathbb M^1\,,\qquad B = \frac{\alpha^2}{2\varepsilon} \mathbb L^1 \bar{DF^{-1}} \bar{\mathbf w}\,,\qquad C = - \frac{1}{2\varepsilon} \bar{DF^{-\top}} \left(\mathbb L^1\right)^\top \,.\]

The accumulation matrix and vector assembled in Accumulator are

\[M = BC \,,\qquad V = B \mathbf V \,.\]

class Variables

Bases: object

property e: FEECVariable

property ions: PICVariable

class Options(solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None)

Bases: object

solver: Literal['pcg', 'cg'] = 'pcg'

precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner'

solver_params: SolverParameters = None

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_coupling.EfieldWeights

Bases: Propagator

Solves the following substep

\[\begin{split}\begin{align} & \frac{\partial \mathbf{E}}{\partial t} = - \frac{\alpha^2}{\varepsilon} \int \mathbf{v} f_1 \, \text{d} \mathbf{v} \,, \\[2mm] & \frac{\partial f_1}{\partial t} = \frac{1}{v_{\text{th}}^2 \varepsilon} \, \mathbf{E} \cdot \mathbf{v} f_0 \,, \end{align}\end{split}\]

which after discretization and in curvilinear coordinates reads

\[\begin{split}\frac{\text{d}}{\text{d} t} w_p &= \frac{f_{0,p}}{s_{0, p}} \frac{1}{v_{\text{th}}^2 \varepsilon} \left[ DF^{-T} (\mathbb{\Lambda}^1)^T \mathbf{e} \right] \cdot \mathbf{v}_p \,, \\[2mm] \frac{\text{d}}{\text{d} t} \mathbb{M}_1 \mathbf{e} &= - \frac{\alpha^2}{\varepsilon} \frac 1N \sum_p w_p \mathbb{\Lambda}^1 \cdot \left( DF^{-1} \mathbf{v}_p \right) \,.\end{split}\]

This is solved using the Crank-Nicolson method

\[\begin{split}\begin{bmatrix} \mathbb{M}_1 \left( \mathbf{e}^{n+1} - \mathbf{e}^n \right) \mathbf{W}^{n+1} - \mathbf{W}^n \end{bmatrix} = \frac{\Delta t}{2} \begin{bmatrix} 0 & - \mathbb{E} \\ \mathbb{W} & 0 \end{bmatrix} \begin{bmatrix} \mathbf{e}^{n+1} + \mathbf{e}^n \\ \mathbf{V}^{n+1} + \mathbf{V}^n \end{bmatrix} \,,\end{split}\]

where

\[\begin{split}\mathbb{E} &= \frac{\alpha^2}{\varepsilon} \frac 1N \mathbb{\Lambda}^1 \cdot \left( DF^{-1} \mathbf{v}_p \right) \,, \\[2mm] \mathbb{W} &= \frac{f_{0,p}}{s_{0,p}} \frac{1}{v_\text{th}^2 \varepsilon} \left( DF^{-1} \mathbf{v}_p \right) \cdot \left(\mathbb{\Lambda}^1\right)^T \,,\end{split}\]

based on the SchurSolver.

The accumulation matrix and vector assembled in Accumulator are

\[BC = \mathbb{E} \mathbb{W} \, , \qquad By_n = \mathbb{E} \mathbf{W} \,.\]

class Variables

Bases: object

property e: FEECVariable

property ions: PICVariable

class Options(alpha: float = 1.0, kappa: float = 1.0, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None)

Bases: object

alpha: float = 1.0

kappa: float = 1.0

solver: Literal['pcg', 'cg'] = 'pcg'

precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner'

solver_params: SolverParameters = None

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_coupling.PressureCoupling6D

Bases: Propagator

FEEC discretization of the following equations: find \(\tilde{\mathbf{U}} \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\) and \(f\) such that

\[\begin{split}\int_\Omega \rho_0 &\frac{\partial \tilde{\mathbf{U}}}{\partial t} \cdot \mathbf V \,\textrm d \mathbf x = - \frac{A_\textnormal{h}}{A_\textnormal{b}} \nabla \cdot \tilde{\mathbb{P}}_{\textnormal{h},\perp} \cdot \mathbf V \,\textrm d \mathbf x \qquad \forall \, \mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\,, \\[2mm] &\frac{\partial f_\textnormal{h}}{\partial t} - \left(\nabla \tilde{\mathbf U}_\perp \cdot \mathbf{v} \right) \cdot \frac{\partial f_\textnormal{h}}{\partial \mathbf{v}} =0\,, \\[2mm] &\tilde{\mathbb{P}}_{\textnormal{h},\perp} = \int_{\mathbb{R}^3}f_\textnormal{h}\mathbf{v}_\perp \mathbf{v}_\perp^\top \,\textnormal{d}^3 \mathbf v\,.\end{split}\]

Time discretization: Crank-Nicolson (implicit mid-point). System size reduction via SchurSolver:

\[\begin{split}\begin{bmatrix} u^{n+1} - u^n \\ V^{n+1} - V^n \end{bmatrix} = \frac{\Delta t}{2} \begin{bmatrix} 0 & (\mathbb M^n)^{-1} V^\top (\bar {\mathcal X})^\top \mathbb G^\top (\bar {\mathbf \Lambda}^1)^\top \bar {DF}^{-1} \\ - {DF}^{-\top} \bar {\mathbf \Lambda}^1 \mathbb G \bar {\mathcal X} V (\mathbb M^n)^{-1} & 0 \end{bmatrix} \begin{bmatrix} {\mathbb M^n}(u^{n+1} + u^n) \\ \bar W (V^{n+1} + V^{n} \end{bmatrix} \,.\end{split}\]

class Variables

Bases: object

property u: FEECVariable

property energetic_ions: PICVariable

class Options(ep_scale: float = 1.0, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None, filter_params: struphy.pic.accumulation.filter.FilterParameters = None, use_perp_model: bool = True)

Bases: object

ep_scale: float = 1.0

u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv'

solver: Literal['pcg', 'cg'] = 'pcg'

precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner'

solver_params: SolverParameters = None

filter_params: FilterParameters = None

use_perp_model: bool = True

allocate()

Allocate all data/objects of the instance.

class GT_MAT_G(derham, MAT, transposed=False)

Bases: LinOpWithTransp

Class for defining LinearOperator corresponding to \(G^\top (\text{MAT}) G \in \mathbb{R}^{3N^0 \times 3N^0}\) where \(\text{MAT} = V^\top (\bar {\mathbf \Lambda}^1)^\top \bar{DF}^{-1} \bar{W} \bar{DF}^{-\top} \bar{\mathbf \Lambda}^1 V \in \mathbb{R}^{3N^1 \times 3N^1}\).

Parameters:

• derham (struphy.feec.psydac_derham.Derham) – Discrete de Rham sequence on the logical unit cube.

• MAT (List of StencilMatrices) – List with six of accumulated pressure terms

property domain

The domain of the linear operator - an element of Vectorspace

property codomain

The codomain of the linear operator - an element of Vectorspace

property dtype

The data type of the coefficients of the linear operator, upon convertion to matrix.

property tosparse

Convert to a sparse matrix in any of the formats supported by scipy.sparse.

property toarray

Convert to Numpy 2D array.

property transposed

transpose()

Transpose the LinearOperator .

If conjugate is True, return the Hermitian transpose.

dot(v, out=None)

dot product between GT_MAT_G and v.

Parameters:

v (StencilVector or BlockVector) – Input FE coefficients from V.coeff_space.

Return type:

A StencilVector or BlockVector from W.coeff_space.

class struphy.propagators.propagators_coupling.CurrentCoupling6DCurrent

Bases: Propagator

FEEC discretization of the following equations: find \(\tilde{\mathbf{U}} \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\) and \(f\) such that

\[\begin{split}\int_\Omega \rho_0 &\frac{\partial \tilde{\mathbf{U}}}{\partial t} \cdot \mathbf V \,\textrm d \mathbf x = - \frac{A_\textnormal{h}}{A_\textnormal{b}} \frac{1}{\varepsilon} \int_\Omega n_\textnormal{h}\mathbf{u}_\textnormal{h} \times(\mathbf{B}_0+\tilde{\mathbf{B}}) \cdot \mathbf V \,\textrm d \mathbf x \qquad \forall \, \mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\,, \\[2mm] &\frac{\partial f_\textnormal{h}}{\partial t} + \frac{1}{\varepsilon} \Big[(\mathbf{B}_0+\tilde{\mathbf{B}})\times\tilde{\mathbf{U}} \Big] \cdot \frac{\partial f_\textnormal{h}}{\partial \mathbf{v}} =0\,, \\[2mm] &n_\textnormal{h}\mathbf{u}_\textnormal{h}=\int_{\mathbb{R}^3}f_\textnormal{h}\mathbf{v}\,\textnormal{d}^3 \mathbf v\,.\end{split}\]

Time discretization: Crank-Nicolson (implicit mid-point). System size reduction via SchurSolver.

class Variables

Bases: object

property ions: PICVariable

property u: FEECVariable

class Options(b_tilde: struphy.models.variables.FEECVariable = None, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None, filter_params: struphy.pic.accumulation.filter.FilterParameters = None, boundary_cut: tuple = (0.0, 0.0, 0.0))

Bases: object

b_tilde: FEECVariable = None

u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv'

solver: Literal['pcg', 'cg'] = 'pcg'

precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner'

solver_params: SolverParameters = None

filter_params: FilterParameters = None

boundary_cut: tuple = (0.0, 0.0, 0.0)

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_coupling.CurrentCoupling5DCurlb

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf U \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\) and \(\mathbf B \in H(\textnormal{div})\) such that

\[\begin{split}\left\{ \begin{aligned} \int \rho_0 &\frac{\partial \tilde{\mathbf U}}{\partial t} \cdot \mathbf V\, \textnormal{d} \mathbf x = - \frac{A_\textnormal{h}}{A_b} \iint \frac{f^\text{vol}}{B^*_\parallel} v_\parallel^2 (\nabla \times \mathbf b_0) \times \mathbf B \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\}\,, \\ &\frac{\partial f}{\partial t} = - \frac{1}{B^*_\parallel} v_\parallel (\nabla \times \mathbf b_0) \cdot (\mathbf B \times \tilde{\mathbf U}) \nabla_{v_\parallel}f \,. \end{aligned} \right.\end{split}\]

Time discretization: Crank-Nicolson (implicit mid-point). System size reduction via SchurSolver:

\[\begin{split}\begin{bmatrix} \mathbf u^{n+1} - \mathbf u^n \\ V_\parallel^{n+1} - V_\parallel^n \end{bmatrix} = \frac{\Delta t}{2} \begin{bmatrix} 0 & - (\mathbb{M}^{2,n})^{-1} \left\{ \mathbb{L}^2 \frac{1}{\bar{\sqrt{g}}} \right\}\cdot_\text{vector} \left\{\bar{b}^{\nabla \times}_0 (\bar{B}^\times_f)^\top \bar{V}_\parallel \frac{1}{\bar{\sqrt{g}}}\right\} \frac{1}{\bar B^{*0}_\parallel}) \\ \frac{1}{\bar B^{*0}_\parallel} \left\{\bar{b}^{\nabla \times}_0 (\bar{B}^\times_f)^\top \bar{V}_\parallel \frac{1}{\bar{\sqrt{g}}}\right\}\, \cdot_\text{vector} \left\{\frac{1}{\bar{\sqrt{g}}}(\mathbb{L}²)^\top\right\} (\mathbb{M}^{2,n})^{-1} & 0 \end{bmatrix} \begin{bmatrix} (\mathbb{M}^{2,n})^{-1} (\mathbf u^{n+1} + \mathbf u^n) \\ \frac{A_\textnormal{h}}{A_b} W (V_\parallel^{n+1} + V_\parallel^n) \end{bmatrix} \,.\end{split}\]

For the detail explanation of the notations, see 2022_DriftKineticCurrentCoupling.

class Variables

Bases: object

property u: FEECVariable

property energetic_ions: PICVariable

class Options(b_tilde: struphy.models.variables.FEECVariable = None, ep_scale: float = 1.0, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None, filter_params: struphy.pic.accumulation.filter.FilterParameters = None)

Bases: object

b_tilde: FEECVariable = None

ep_scale: float = 1.0

u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv'

solver: Literal['pcg', 'cg'] = 'pcg'

precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner'

solver_params: SolverParameters = None

filter_params: FilterParameters = None

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_coupling.CurrentCoupling5DGradB

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf U \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\) and \(\mathbf B \in H(\textnormal{div})\) such that

\[\begin{split}\left\{ \begin{aligned} \int \rho_0 &\frac{\partial \tilde{\mathbf U}}{\partial t} \cdot \mathbf V \, \textnormal{d} \mathbf x = - \frac{A_\textnormal{h}}{A_b} \iint \mu \frac{f^\text{vol}}{B^*_\parallel} (\mathbf b_0 \times \nabla B_\parallel) \times \mathbf B \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\}\,, \\ &\frac{\partial f}{\partial t} = \frac{1}{B^*_\parallel} \left[ \mathbf b_0 \times (\tilde{\mathbf U} \times \mathbf B) \right] \cdot \nabla f\,. \end{aligned} \right.\end{split}\]

Time discretization: Explicit (‘rk4’, ‘forward_euler’, ‘heun2’, ‘rk2’, ‘heun3’).

\[\begin{split}\begin{bmatrix} \dot{\mathbf u}\\ \dot{\mathbf H} \end{bmatrix} = \begin{bmatrix} 0 & (\mathbb{M}^{2,n})^{-1} \mathbb{L}² \frac{1}{\bar{\sqrt{g}}} \frac{1}{\bar B^{*0}_\parallel}\bar{B}^\times_f \bar{G}^{-1} \bar{b}^\times_0 \bar{G}^{-1} \\ -\bar{G}^{-1} \bar{b}^\times_0 \bar{G}^{-1} \bar{B}^\times_f \frac{1}{\bar B^{*0}_\parallel} \frac{1}{\bar{\sqrt{g}}} (\mathbb{L}²)^\top (\mathbb{M}^{2,n})^{-1} & 0 \end{bmatrix} \begin{bmatrix} \mathbb M^{2,n} \mathbf u \\ \frac{A_\textnormal{h}}{A_b} \bar M \bar W \overline{\nabla B}_\parallel \end{bmatrix} \,.\end{split}\]

For the detail explanation of the notations, see 2022_DriftKineticCurrentCoupling.

class Variables

Bases: object

property u: FEECVariable

property energetic_ions: PICVariable

class Options(b_tilde: struphy.models.variables.FEECVariable = None, ep_scale: float = 1.0, algo: Literal['discrete_gradient', 'explicit'] = 'explicit', butcher: struphy.ode.utils.ButcherTableau = None, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None, filter_params: struphy.pic.accumulation.filter.FilterParameters = None, dg_solver_params: struphy.linear_algebra.solver.DiscreteGradientSolverParameters = None)

Bases: object

OptsAlgo

alias of Literal[‘discrete_gradient’, ‘explicit’]

b_tilde: FEECVariable = None

ep_scale: float = 1.0

algo: Literal['discrete_gradient', 'explicit'] = 'explicit'

butcher: ButcherTableau = None

u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv'

solver: Literal['pcg', 'cg'] = 'pcg'

precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner'

solver_params: SolverParameters = None

filter_params: FilterParameters = None

dg_solver_params: DiscreteGradientSolverParameters = None

allocate()

Allocate all data/objects of the instance.