Field propagators

Only FEEC variables are updated.

class struphy.propagators.propagators_fields.Maxwell

Bases: Propagator

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

\[\begin{split}&\int_\Omega \frac{\partial \mathbf E}{\partial t} \cdot \mathbf F \, \textrm d \mathbf x - \int_\Omega \mathbf B \cdot \nabla \times \mathbf F \,\textrm d \mathbf x = 0\,, \qquad \forall \, \mathbf F \in H(\textnormal{curl}) \,. \\[2mm] &\frac{\partial \mathbf B}{\partial t} + \nabla\times\mathbf E = 0\,.\end{split}\]

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

class Variables

Bases: object

property e: FEECVariable

property b: FEECVariable

class Options(algo: Literal['implicit', 'explicit'] = 'implicit', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None, butcher: struphy.ode.utils.ButcherTableau = None)

Bases: object

OptsAlgo

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

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

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

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

solver_params: SolverParameters = None

butcher: ButcherTableau = None

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_fields.OhmCold

Bases: Propagator

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

\[\begin{split}\int_\Omega \frac{1}{n_0} \frac{\partial \mathbf j}{\partial t} \cdot \mathbf F \,\textrm d \mathbf x &= \frac{1}{\varepsilon} \int_\Omega \mathbf E \cdot \mathbf F \,\textrm d \mathbf x \qquad \forall \,\mathbf F \in H(\textnormal{curl})\,, \\[2mm] -\frac{\partial \mathbf E}{\partial t} &= \frac{\alpha^2}{\varepsilon} \mathbf j \,,\end{split}\]

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

\[\begin{split}\begin{bmatrix} \mathbf j^{n+1} - \mathbf j^n \\ \mathbf e^{n+1} - \mathbf e^n \end{bmatrix} = \frac{\Delta t}{2} \begin{bmatrix} 0 & \frac{1}{\varepsilon} \mathbb M_{1/n_0}^{-1} \\ - \frac{1}{\varepsilon} \mathbb M_{1/n_0}^{-1} & 0 \end{bmatrix} \begin{bmatrix} \alpha^2 \mathbb M_{1/n_0} (\mathbf j^{n+1} + \mathbf j^{n}) \\ \mathbb M_1 (\mathbf e^{n+1} + \mathbf e^{n}) \end{bmatrix} \,.\end{split}\]

class Variables

Bases: object

property j: FEECVariable

property e: FEECVariable

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_fields.JxBCold

Bases: Propagator

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

\[\int_\Omega \frac{1}{n_0} \frac{\partial \mathbf j}{\partial t} \cdot \mathbf F \,\textrm d \mathbf x = \frac{1}{\varepsilon} \int_\Omega \frac{1}{n_0} (\mathbf j \times \mathbf B_0) \cdot \mathbf F \,\textrm d \mathbf x \qquad \forall \,\mathbf F \in H(\textnormal{curl})\,,\]

Time discretization: Crank-Nicolson (implicit mid-point), such that

\[\mathbb M_{1/n_0} \left( \mathbf j^{n+1} - \mathbf j^n \right) = \frac{\Delta t}{2} \frac{1}{\varepsilon} \mathbb M_{B_0/n_0} \left( \mathbf j^{n+1} - \mathbf j^n \right)\,.\]

class Variables

Bases: object

property j: FEECVariable

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_fields.ShearAlfven

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{align}\begin{aligned}&\int_\Omega \rho_0\frac{\partial \tilde{\mathbf{U}}}{\partial t} \cdot \mathbf V\,\textnormal d \mathbf x =\int_\Omega \tilde{\mathbf B } \cdot \nabla \times (\mathbf{B}_0 \times \mathbf V) \,\textnormal d \mathbf x \qquad \forall \,\mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\,,\\&\frac{\partial \tilde{\mathbf{B}}}{\partial t} - \nabla\times(\tilde{\mathbf{U}} \times \mathbf{B}_0) = 0\,.\end{aligned}\end{align} \]

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

\[\begin{split}\begin{bmatrix} \mathbf u^{n+1} - \mathbf u^n \\ \mathbf b^{n+1} - \mathbf b^n \end{bmatrix} = \frac{\Delta t}{2} \begin{bmatrix} 0 & (\mathbb M^\rho_\alpha)^{-1} \mathcal {T^\alpha}^\top \mathbb C^\top \\ - \mathbb C \mathcal {T^\alpha} (\mathbb M^\rho_\alpha)^{-1} & 0 \end{bmatrix} \begin{bmatrix} {\mathbb M^\rho_\alpha}(\mathbf u^{n+1} + \mathbf u^n) \\ \mathbb M_2(\mathbf b^{n+1} + \mathbf b^n) \end{bmatrix} ,\end{split}\]

where \(\alpha \in \{1, 2, v\}\) and \(\mathbb M^\rho_\alpha\) is a weighted mass matrix in \(\alpha\)-space, the weight being \(\rho_0\), the MHD equilibirum density. The solution of the above system is based on the Schur complement.

class Variables

Bases: object

property u: FEECVariable

property b: FEECVariable

class Options(u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', algo: Literal['implicit', 'explicit'] = 'implicit', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None, butcher: struphy.ode.utils.ButcherTableau = None)

Bases: object

OptsAlgo

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

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

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

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

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

solver_params: SolverParameters = None

butcher: ButcherTableau = None

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_fields.ShearAlfvenB1

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{curl})\) such that

\[\begin{split}&\int_\Omega \rho_0\frac{\partial \tilde{\mathbf{U}}}{\partial t} \cdot \mathbf V\,\textnormal d \mathbf x =\int_\Omega \tilde{\mathbf B } \cdot \nabla \times (\mathbf{B}_0 \times \mathbf V) \,\textnormal d \mathbf x \qquad \forall \,\mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\,, \\[2mm] &\int_\Omega \frac{\partial \tilde{\mathbf{B}}}{\partial t} \cdot \mathbf C\,\textnormal d \mathbf x - \int_\Omega \mathbf C \cdot\nabla\times \left( \tilde{\mathbf{U}} \times \mathbf{B}_0 \right) \textrm d \mathbf x = 0 \qquad \forall \, \mathbf C \in H(\textrm{curl})\,.\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 \\ \mathbf b^{n+1} - \mathbf b^n \end{bmatrix} = \frac{\Delta t}{2} \begin{bmatrix} 0 & (\mathbb M^\rho_2)^{-1} \mathcal {T^2}^\top \mathbb C \mathbb M_1^{-1}\\ - \mathbb M_1^{-1} \mathbb C^\top \mathcal {T^2} (\mathbb M^\rho_2)^{-1} & 0 \end{bmatrix} \begin{bmatrix} {\mathbb M^\rho_2}(\mathbf u^{n+1} + \mathbf u^n) \\ \mathbb M_1(\mathbf b^{n+1} + \mathbf b^n) \end{bmatrix} ,\end{split}\]

where \(\mathbb M^\rho_2\) is a weighted mass matrix in 2-space, the weight being \(\rho_0\), the MHD equilibirum density.

class Variables

Bases: object

property u: FEECVariable

property b: FEECVariable

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

Bases: object

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

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

solver_params: SolverParameters = None

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

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

solver_params_M1: SolverParameters = None

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_fields.Hall

Bases: Propagator

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

\[\int_\Omega \frac{\partial \tilde{\mathbf{B}}}{\partial t} \cdot \mathbf C\,\textnormal d \mathbf x + \frac{1}{\varepsilon} \int_\Omega \nabla\times \mathbf C \cdot \left( \frac{\nabla\times \tilde{\mathbf{B}}}{\rho_0}\times \mathbf{B}_0 \right) \textrm d \mathbf x = 0 \qquad \forall \, \mathbf C \in H(\textrm{curl})\]

Time discretization: Crank-Nicolson (implicit mid-point):

\[\mathbf b^{n+1} - \mathbf b^n = \frac{\Delta t}{2} \mathbb M_1^{-1} \mathbb C^\top \mathbb M^{\mathcal{T},\rho}_2 \mathbb C (\mathbf b^{n+1} + \mathbf b^n) ,\]

where \(\mathbb M^{\mathcal{T},\rho}_2\) is a weighted mass matrix in 2-space, the weight being \(\frac{\mathcal{T}}{\rho_0}\), the MHD equilibirum density \(\rho_0\) as a 0-form, and rotation matrix \(\mathcal{T} \vec v = \vec B^2_{\textnormal{eq}} \times \vec v\,,\). The solution of the above system is based on the Pre-conditioned Biconjugate Gradient Stabilized algortihm (PBiConjugateGradientStab).

class Variables

Bases: object

property b: FEECVariable

class Options(solver: Literal['pbicgstab', 'bicgstab', 'GMRES'] = 'pbicgstab', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None, epsilon_from: struphy.models.species.Species = None)

Bases: object

solver: Literal['pbicgstab', 'bicgstab', 'GMRES'] = 'pbicgstab'

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

solver_params: SolverParameters = None

epsilon_from: Species = None

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_fields.Magnetosonic

Bases: Propagator

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

\[ \begin{align}\begin{aligned}&\frac{\partial \tilde \rho}{\partial t}+\nabla\cdot(\rho_0 \tilde{\mathbf{U}})=0\,,\\\int \rho_0&\frac{\partial \tilde{\mathbf{U}}}{\partial t} \cdot \mathbf V\,\textrm d \mathbf x - \int \tilde p\, \nabla \cdot \mathbf V \,\textrm d \mathbf x =\int (\nabla\times\mathbf{B}_0)\times \tilde{\mathbf{B}} \cdot \mathbf V\,\textrm d \mathbf x \qquad \forall \ \mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\,,\\&\frac{\partial \tilde p}{\partial t} + \nabla\cdot(p_0 \tilde{\mathbf{U}}) + \frac{2}{3}\,p_0\nabla\cdot \tilde{\mathbf{U}}=0\,.\end{aligned}\end{align} \]

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

\[\begin{split}\begin{bmatrix} \mathbf u^{n+1} - \mathbf u^n \\ \mathbf p^{n+1} - \mathbf p^n \end{bmatrix} = \frac{\Delta t}{2} \begin{bmatrix} 0 & (\mathbb M^\rho_\alpha)^{-1} {\mathcal U^\alpha}^\top \mathbb D^\top \mathbb M_3 \\ - \mathbb D \mathcal S^\alpha - (\gamma - 1) \mathcal K^\alpha \mathbb D \mathcal U^\alpha & 0 \end{bmatrix} \begin{bmatrix} (\mathbf u^{n+1} + \mathbf u^n) \\ (\mathbf p^{n+1} + \mathbf p^n) \end{bmatrix} + \begin{bmatrix} \Delta t (\mathbb M^\rho_\alpha)^{-1} \mathbb M^J_\alpha \mathbf b^n \\ 0 \end{bmatrix},\end{split}\]

where \(\alpha \in \{1, 2, v\}\) and \(\mathcal U^2 = \mathbb Id\); moreover, \(\mathbb M^\rho_\alpha\) and \(\mathbb M^J_\alpha\) are weighted mass matrices in \(\alpha\)-space, the weights being the MHD equilibirum density \(\rho_0\) and the curl of the MHD equilibrium current density \(\mathbf J_0 = \nabla \times \mathbf B_0\). Density update is decoupled:

\[\boldsymbol{\rho}^{n+1} = \boldsymbol{\rho}^n - \frac{\Delta t}{2} \mathbb D \mathcal Q^\alpha (\mathbf u^{n+1} + \mathbf u^n) \,.\]

class Variables

Bases: object

property n: FEECVariable

property u: FEECVariable

property p: FEECVariable

class Options(b_field: struphy.models.variables.FEECVariable = None, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', solver: Literal['pbicgstab', 'bicgstab', 'GMRES'] = 'pbicgstab', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None)

Bases: object

b_field: FEECVariable = None

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

solver: Literal['pbicgstab', 'bicgstab', 'GMRES'] = 'pbicgstab'

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

solver_params: SolverParameters = None

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_fields.MagnetosonicUniform

Bases: Propagator

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

\[ \begin{align}\begin{aligned}&\frac{\partial \tilde \rho}{\partial t}+\nabla\cdot(\rho_0 \tilde{\mathbf{U}})=0\,,\\\int \rho_0&\frac{\partial \tilde{\mathbf{U}}}{\partial t} \cdot \mathbf V\,\textrm d \mathbf x - \int \tilde p\, \nabla \cdot \mathbf V \,\textrm d \mathbf x = 0 \qquad \forall \ \mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\,,\\&\frac{\partial \tilde p}{\partial t} + \frac{5}{3}\,p_0\nabla\cdot \tilde{\mathbf{U}}=0\,.\end{aligned}\end{align} \]

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

\[\begin{split}\begin{bmatrix} \mathbf u^{n+1} - \mathbf u^n \\ \mathbf p^{n+1}_i - \mathbf p^n_i \end{bmatrix} = \frac{\Delta t}{2} \begin{bmatrix} 0 & (\mathbb M^\rho_2)^{-1} \mathbb D^\top \mathbb M_3 \\ - \gamma \mathcal K^3 \mathbb D & 0 \end{bmatrix} \begin{bmatrix} (\mathbf u^{n+1} + \mathbf u^n) \\ (\mathbf p^{n+1}_i + \mathbf p^n_i) \end{bmatrix} ,\end{split}\]

where \(\mathbb M^\rho_2\) is a weighted mass matrix in 2-space, the weight being the MHD equilibirum density \(\rho_0\). Furthermore, \(\mathcal K^3\) is the basis projection operator given by :

\[\mathcal{K}^3_{ijk,mno} := \hat{\Pi}^3_{ijk} \left[ \frac{\hat{p}^3_{\text{eq}}}{\sqrt{g}}\Lambda^3_{mno} \right] \,.\]

The solution of the above system is based on the Schur complement.

Decoupled density update:

\[\boldsymbol{\rho}^{n+1} = \boldsymbol{\rho}^n - \frac{\Delta t}{2} \mathcal Q \mathbb D (\mathbf u^{n+1} + \mathbf u^n) \,.\]

Parameters:

• n (psydac.linalg.stencil.StencilVector) – FE coefficients of a discrete 3-form.

• u (psydac.linalg.block.BlockVector) – FE coefficients of MHD velocity 2-form.

• p (psydac.linalg.stencil.StencilVector) –

  FE coefficients of a discrete 3-form.

  **paramsdict

  Solver- and/or other parameters for this splitting step.

class Variables

Bases: object

property n: FEECVariable

property u: FEECVariable

property p: FEECVariable

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

Bases: object

solver: Literal['pbicgstab', 'bicgstab', 'GMRES'] = 'pbicgstab'

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

solver_params: SolverParameters = None

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_fields.FaradayExtended(a, **params)

Bases: Propagator

Equations: Faraday’s law

\[\begin{split}\begin{align*} & \frac{\partial {\mathbf A}}{\partial t} = - \frac{\nabla \times (\nabla \times {\mathbf A} + {\mathbf B}_0) }{n} \times (\nabla \times {\mathbf A} + {\mathbf B}_0) - \frac{\int ({\mathbf A} - {\mathbf p}f \mathrm{d}{\mathbf p})}{n} \times (\nabla \times {\mathbf A} + {\mathbf B}_0), \\ & n = \int f \mathrm{d}{\mathbf p}. \end{align*}\end{split}\]

Mid-point rule:

\[\begin{split}\begin{align*} & \left[ \mathbb{M}_1 - \frac{\Delta t}{2} \mathbb{F}(\hat{n}^0_h, {\mathbf a}^{n+\frac{1}{2}}) \mathbb{M}_1^{-1} (\mathbb{P}_1^\top \mathbb{W} \mathbb{P}_1 + \mathbb{C}^\top \mathbb{M}_2 \mathbb{C} ) \right] {\mathbf a}^{n+1} \\ & = \mathbb{M}_1 {\mathbf a}^n + \frac{\Delta t}{2} \mathbb{F}(\hat{n}^0_h, {\mathbf a}^{n+\frac{1}{2}}) \mathbb{M}_1^{-1} (\mathbb{P}_1^\top \mathbb{W} \mathbb{P}_1 + \mathbb{C}^\top \mathbb{M}_2 \mathbb{C} ) {\mathbf a}^{n+1} \\ & - \Delta t \mathbb{F}(\hat{n}^0_h, {\mathbf a}^{n+\frac{1}{2}}) \mathbb{M}_1^{-1} \mathbb{P}_1^\top \mathbb{W} {\mathbf P}^n\\ & + \Delta t \mathbb{F}(\hat{n}^0_h, {\mathbf a}^{n+\frac{1}{2}}) \mathbb{M}_1^{-1} \mathbb{C}^\top \mathbb{M}_2 {\mathbf b}_0\\ & \mathbb{F}_{ij} = - \int \frac{1}{\hat{n}^0_h \sqrt{g}} G (\nabla \times {\mathbf A} + {\mathbf B}_0) \cdot (\Lambda^1_i \times \Lambda^1_j) \mathrm{d}{\boldsymbol \eta}. \end{align*}\end{split}\]

Parameters:

• a (psydac.linalg.block.BlockVector) – FE coefficients of vector potential.

• **params (dict) – Solver- and/or other parameters for this splitting step.

class struphy.propagators.propagators_fields.CurrentCoupling6DDensity

Bases: Propagator

FEEC discretization of the following equations: find \(\tilde{\mathbf{U}} \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\) 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}\tilde{\mathbf{U}} \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] &n_\textnormal{h}=\int_{\mathbb{R}^3}f_\textnormal{h}\,\textnormal{d}^3 \mathbf v\,.\end{split}\]

Time discretization: Crank-Nicolson (implicit mid-point).

class Variables

Bases: object

property u: FEECVariable

class Options(energetic_ions: struphy.models.variables.PICVariable = None, 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

energetic_ions: PICVariable = None

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_fields.ShearAlfvenCurrentCoupling5D

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} = \int \left(\tilde{\mathbf B} - \frac{A_\textnormal{h}}{A_b} \iint f^\text{vol} \mu \mathbf{b}_0\textnormal{d} \mu \textnormal{d} v_\parallel \right) \cdot \nabla \times (\tilde{\mathbf B} \times \mathbf V) \, \textnormal{d} \mathbf{x} \quad \forall \, \mathbf V \in \{H(\textnormal{curl}), H(\textnormal{div}), (H^1)^3\}\,, \,, \\ &\frac{\partial \tilde{\mathbf B}}{\partial t} = - \nabla \times (\tilde{\mathbf B} \times \tilde{\mathbf U}) \,. \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 \\ \mathbf b^{n+1} - \mathbf b^n \end{bmatrix} = \frac{\Delta t}{2} \,. \begin{bmatrix} 0 & (\mathbb M^{2,n})^{-1} \mathcal {T^2}^\top \mathbb C^\top \\ - \mathbb C \mathcal {T^2} (\mathbb M^{2,n})^{-1} & 0 \end{bmatrix} \begin{bmatrix} {\mathbb M^{2,n}}(\mathbf u^{n+1} + \mathbf u^n) \\ \mathbb M_2(\mathbf b^{n+1} + \mathbf b^n) + \sum_k^{N_p} \omega_k \mu_k \hat{\mathbf b}¹_0 (\boldsymbol \eta_k) \cdot \left(\frac{1}{\sqrt{g(\boldsymbol \eta_k)}} \vec \Lambda² (\boldsymbol \eta_k) \right) \end{bmatrix} \,,\end{split}\]

where \(\mathcal{T}^2 = \hat \Pi \left[\frac{\tilde{\mathbf B}^2}{\sqrt{g} \times \vec \Lambda^2\right]\) and \(\mathbb M^{2,n}\) is a WeightedMassOperators being weighted with \(\rho_\text{eq}\), the MHD equilibirum density.

class Variables

Bases: object

property u: FEECVariable

property b: FEECVariable

class Options(energetic_ions: struphy.models.variables.PICVariable = None, ep_scale: float = 1.0, u_space: Literal['Hcurl', 'Hdiv', 'H1vec'] = 'Hdiv', algo: Literal['implicit', 'explicit'] = 'implicit', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixDiagonalPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None, filter_params: struphy.pic.accumulation.filter.FilterParameters = None, butcher: struphy.ode.utils.ButcherTableau = None, nonlinear: bool = True)

Bases: object

OptsAlgo

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

energetic_ions: PICVariable = None

ep_scale: float = 1.0

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

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

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

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

solver_params: SolverParameters = None

filter_params: FilterParameters = None

butcher: ButcherTableau = None

nonlinear: bool = True

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_fields.CurrentCoupling5DDensity

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

\[\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:

\[\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.

class Variables

Bases: object

property u: FEECVariable

class Options(energetic_ions: struphy.models.variables.PICVariable = None, 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

energetic_ions: PICVariable = None

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_fields.ImplicitDiffusion

Bases: Propagator

Weak, implicit discretization of the diffusion (or heat) equation (can be used as a Poisson solver too).

Find \(\phi \in H^1\) such that

\[\int_\Omega \psi\, n_0(\mathbf x)\frac{\partial \phi}{\partial t}\,\textrm d \mathbf x + \int_\Omega \nabla \psi^\top D_0(\mathbf x) \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 \(n_0, \rho_i:\Omega \to \mathbb R\) are real-valued functions and \(D_0:\Omega \to \mathbb R^{3\times 3}\) is a positive diffusion matrix. Boundary terms from integration by parts are assumed to vanish. The equation is discretized as

\[\left( \frac{\sigma_1}{\Delta t} \mathbb M^0_{n_0} + \mathbb G^\top \mathbb M^1_{D_0} \mathbb G \right)\, \boldsymbol\phi^{n+1} = \frac{\sigma_2}{\Delta t} \mathbb M^0_{n_0} \boldsymbol\phi^{n} + \frac{\sigma_3}{\Delta t} \sum_i(\Lambda^0, \rho_i )_{L^2}\,,\]

where \(M^0_{n_0}\) and \(M^1_{D_0}\) are WeightedMassOperators and \(\sigma_1, \sigma_2, \sigma_3 \in \mathbb R\) are artificial parameters that can be tuned to change the model (see Notes).

Notes

• \(\sigma_1=\sigma_2=0\) and \(\sigma_3 = \Delta t\): Poisson solver with a given charge density \(\sum_i\rho_i\).

• \(\sigma_2=0\) and \(\sigma_1 = \sigma_3 = \Delta t\) : Poisson with adiabatic electrons.

• \(\sigma_1=\sigma_2=1\) and \(\sigma_3 = 0\): Implicit heat equation.

class Variables

Bases: object

property phi: FEECVariable

class Options(sigma_1: float = 1.0, sigma_2: float = 0.0, sigma_3: float = 1.0, divide_by_dt: bool = False, stab_mat: Literal['M0', 'M0ad', 'Id'] = 'M0', diffusion_mat: Literal['M1', 'M1perp'] = 'M1', rho: struphy.models.variables.FEECVariable | Callable | tuple[struphy.pic.accumulation.particles_to_grid.AccumulatorVector, struphy.pic.base.Particles] | list = None, rho_coeffs: float | list = None, x0: psydac.linalg.stencil.StencilVector = None, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None)

Bases: object

OptsStabMat

alias of Literal[‘M0’, ‘M0ad’, ‘Id’]

OptsDiffusionMat

alias of Literal[‘M1’, ‘M1perp’]

sigma_1: float = 1.0

sigma_2: float = 0.0

sigma_3: float = 1.0

divide_by_dt: bool = False

stab_mat: Literal['M0', 'M0ad', 'Id'] = 'M0'

diffusion_mat: Literal['M1', 'M1perp'] = 'M1'

rho: FEECVariable | Callable | tuple[AccumulatorVector, Particles] | list = None

rho_coeffs: float | list = None

x0: StencilVector = None

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

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

solver_params: SolverParameters = None

allocate()

Allocate all data/objects of the instance.

property sources: list[StencilVector | FEECVariable | AccumulatorVector]

Right-hand side of the equation (sources).

property coeffs: list[float]

Same length as self.sources. Coefficients multiplied with sources before solve (default is 1.0).

property x0

psydac.linalg.stencil.StencilVector or struphy.polar.basic.PolarVector. First guess of the iterative solver.

class struphy.propagators.propagators_fields.Poisson

Bases: ImplicitDiffusion

Weak discretization of the (stabilized) Poisson equation.

Find \(\phi \in H^1\) such that

\[\epsilon \int_\Omega \psi\, \phi\,\textrm d \mathbf x + \int_\Omega \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 \(\epsilon \in \mathbb R\) is a stabilization parameter. Boundary terms from integration by parts are assumed to vanish. The equation is discretized as

\[\left( \epsilon\,\mathbb S + \mathbb G^\top \mathbb M^1 \mathbb G \right)\, \boldsymbol\phi^{n+1} = \sum_i(\Lambda^0, \rho_i )_{L^2}\,,\]

where \(\mathbb M^1\) is the \(H(\textnormal{curl})\)-mass matrix and \(\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 AccumulatorVector, and the second entry must be 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).

class Options(stab_eps: float = 0.0, stab_mat: Literal['M0', 'M0ad', 'Id'] = 'Id', rho: struphy.models.variables.FEECVariable | Callable | tuple[struphy.pic.accumulation.particles_to_grid.AccumulatorVector, struphy.pic.base.Particles] | list = None, rho_coeffs: float | list = None, x0: psydac.linalg.stencil.StencilVector = None, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None)

Bases: object

OptsStabMat

alias of Literal[‘M0’, ‘M0ad’, ‘Id’]

stab_eps: float = 0.0

stab_mat: Literal['M0', 'M0ad', 'Id'] = 'Id'

rho: FEECVariable | Callable | tuple[AccumulatorVector, Particles] | list = None

rho_coeffs: float | list = None

x0: StencilVector = None

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

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

solver_params: SolverParameters = None

class struphy.propagators.propagators_fields.VariationalMomentumAdvection

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf u \in (H^1)^3\) such that

\[\int_{\Omega} \partial_t ( \rho \mathbf{u}) \cdot \mathbf{v} \,\textrm d \mathbf x - \int_{\Omega} \rho \mathbf{u} \cdot [\mathbf{u}, \mathbf{v}] \, \textrm d \mathbf x = 0 \,.\]

On the logical domain:

\[\int_{\hat{\Omega}} \partial_t ( \hat{\rho}^3 \hat{\mathbf{u}}) \cdot G \hat{\mathbf{v}} \,\textrm d \boldsymbol \eta - \int_{\hat{\Omega}} \hat{\rho}^3 \hat{\mathbf{u}} \cdot G [\hat{\mathbf{u}}, \hat{\mathbf{v}}] \, \textrm d \boldsymbol \eta = 0 \,,\]

which is discretized as

\[\mathbb M^v[\hat{\rho}_h^n] \frac{\mathbf u^{n+1}- \mathbf u^n}{\Delta t} - \left(\sum_{\mu} (\hat{\Pi}^{0}[\hat{\mathbf u}_h^{n+1/2} \cdot \vec{\boldsymbol \Lambda}^1] \mathbb G P_{\mu} - \hat{\Pi}^0[\hat{\mathbf A}^1_{\mu,h} \cdot \vec{\boldsymbol \Lambda}^v])^\top P_i \right) \mathbb M^v[\hat{\rho}_h^n] \mathbf u^{n} = 0 ~ .\]

where \(P_\mu\) stand for the CoordinateProjector and the weights in the the two BasisProjectionOperator and the WeightedMassOperator are given by

\[\hat{\mathbf{u}}_h^{n+1/2} = (\mathbf{u}^{n+1/2})^\top \vec{\boldsymbol \Lambda}^v \in (V_h^0)^3 \,, \qquad \hat{\mathbf A}^1_{\mu,h} = \nabla P_\mu((\mathbf u^{n+1/2})^\top \vec{\boldsymbol \Lambda}^v)] \in V_h^1\,, \qquad \hat{\rho}_h^{n} = (\rho^{n})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \,.\]

class Variables

Bases: object

property u: FEECVariable

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

Bases: object

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

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

solver_params: SolverParameters = None

nonlin_solver: NonlinearSolverParameters = None

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_fields.VariationalDensityEvolve

Bases: Propagator

FEEC discretization of the following equations: find \(\rho \in L^2\) and \(\mathbf u \in (H^1)^3\) such that

\[\begin{split}&\partial_t \rho + \nabla \cdot ( \tilde{\rho} \mathbf u ) = 0 \,, \\[4mm] &\int_\Omega \partial_t (\rho \mathbf u) \cdot \mathbf v\,\textrm d \mathbf x + \int_\Omega \left(\frac{|\mathbf u|^2}{2} - \frac{\partial(\rho \mathcal U(\rho))}{\partial \rho}\right) \nabla \cdot (\tilde{\rho} \mathbf v) \,\textrm d \mathbf x = 0 \qquad \forall \, \mathbf v \in (H^1)^3\,.\end{split}\]

Where \(\tilde{\rho}\) is either \(\rho\) for full-f models, \(\rho_0\) for linear models or \(\rho_0+\rho\) for \(\delta f\) models.

In the case of linear model, the second equation is not updated.

On the logical domain:

\[\begin{split}\begin{align} &\partial_t \hat{\rho}^3 + \nabla \cdot ( \hat{\rho}^3 \hat{\mathbf{u}} ) = 0 ~ , \\[4mm] &\int_{\hat{\Omega}} \partial_t ( \hat{\rho}^3 \hat{\mathbf{u}}) \cdot G \hat{\mathbf{v}} \, \textrm d \boldsymbol \eta + \int_{\hat{\Omega}} \left( \frac{| DF \hat{\mathbf{u}} |^2}{2} - \frac{\partial (\hat{\rho}^3 \mathcal U)}{\partial \hat{\rho}^3} \right) \nabla \cdot (\hat{\rho}^3 \hat{\mathbf{v}}) \, \textrm d \boldsymbol \eta = 0 ~ , \\[2mm] \end{align}\end{split}\]

where \(\mathcal U\) depends on the chosen model. It is discretized as

\[\begin{split}\begin{align} &\frac{\mathbb M^v[\hat{\rho}_h^{n+1}] \mathbf u^{n+1}- \mathbb M^v[\hat{\rho}_h^n] \mathbf u^n}{\Delta t} + (\mathbb D \hat{\Pi}^{2}[\hat{\rho_h^{n}} \vec{\boldsymbol \Lambda}^v])^\top \hat{l}^3\left(\frac{DF \hat{\mathbf{u}}_h^{n+1} \cdot DF \hat{\mathbf{u}}_h^{n}}{2} - \frac{\hat{\rho}_h^{n+1}\mathcal U(\hat{\rho}_h^{n+1})-\hat{\rho}_h^{n}\mathcal U(\hat{\rho}_h^{n})}{\hat{\rho}_h^{n+1}-\hat{\rho}_h^n} \right) = 0 ~ , \\[2mm] &\frac{\boldsymbol \rho^{n+1}- \boldsymbol \rho^n}{\Delta t} + \mathbb D \hat{\Pi}^{2}[\hat{\rho_h^{n}} \vec{\boldsymbol \Lambda}^v] \mathbf u^{n+1/2} = 0 ~ , \\[2mm] \end{align}\end{split}\]

where \(\hat{l}^3(f)\) denotes the vector representing the linear form \(v_h \mapsto \int_{\hat{\Omega}} f(\boldsymbol \eta) v_h(\boldsymbol \eta) d \boldsymbol \eta\), that is the vector with components

\[\hat{l}^3(f)_{ijk}=\int_{\hat{\Omega}} f \Lambda^3_{ijk} \textrm d \boldsymbol \eta\]

and the weights in the the BasisProjectionOperator and the WeightedMassOperator are given by

\[\hat{\mathbf{u}}_h^{k} = (\mathbf{u}^{k})^\top \vec{\boldsymbol \Lambda}^v \in (V_h^0)^3 \, \text{for k in} \{n, n+1/2, n+1\}, \qquad \hat{\rho}_h^{k} = (\rho^{k})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \, \text{for k in} \{n, n+1/2, n+1\} .\]

class Variables

Bases: object

property rho: FEECVariable

property u: FEECVariable

class Options(model: Literal['pressureless', 'barotropic', 'full', 'full_p', 'full_q', 'linear', 'deltaf', 'linear_q', 'deltaf_q'] = 'barotropic', gamma: float = 1.6666666666666667, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None, nonlin_solver: struphy.linear_algebra.solver.NonlinearSolverParameters = None, s: struphy.models.variables.FEECVariable = None)

Bases: object

OptsModel

alias of Literal[‘pressureless’, ‘barotropic’, ‘full’, ‘full_p’, ‘full_q’, ‘linear’, ‘deltaf’, ‘linear_q’, ‘deltaf_q’]

model: Literal['pressureless', 'barotropic', 'full', 'full_p', 'full_q', 'linear', 'deltaf', 'linear_q', 'deltaf_q'] = 'barotropic'

gamma: float = 1.6666666666666667

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

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

solver_params: SolverParameters = None

nonlin_solver: NonlinearSolverParameters = None

s: FEECVariable = None

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_fields.VariationalEntropyEvolve

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf u \in (H^1)^3\) and \(s \in L^2\) such that

\[\begin{split}&\int_\Omega \partial_t (\rho \mathbf u) \cdot \mathbf v\,\textrm d \mathbf x - \int_\Omega \frac{\partial(\rho \mathcal U)}{\partial s} \nabla \cdot (s \mathbf v) \,\textrm d \mathbf x = 0 \qquad \forall \, \mathbf v \in (H^1)^3\,, \\[4mm] &\partial_t s + \nabla \cdot ( s \mathbf u ) = 0 \,.\end{split}\]

On the logical domain:

\[\begin{split}\begin{align} &\int_{\hat{\Omega}} \partial_t ( \hat{\rho}^3 \hat{\mathbf{u}}) \cdot G \hat{\mathbf{v}} \, \textrm d \boldsymbol \eta - \int_{\hat{\Omega}} \left(\frac{\partial \hat{\rho}^3 \mathcal U}{\partial \hat{s}} \right) \nabla \cdot (\hat{s} \hat{\mathbf{v}}) \, \textrm d \boldsymbol \eta = 0 ~ , \\[2mm] &\partial_t \hat{s} + \nabla \cdot ( \hat{s} \hat{\mathbf{u}} ) = 0 ~ , \end{align}\end{split}\]

where \(\mathcal U\) depends on the chosen model. It is discretized as

\[\begin{split}\begin{align} &\mathbb M^v[\hat{\rho}_h^{n}] \frac{ \mathbf u^{n+1}-\mathbf u^n}{\Delta t} - (\mathbb D \hat{\Pi}^{2}[\hat{s_h^{n}} \vec{\boldsymbol \Lambda}^v])^\top \hat{l}^3\left( \frac{\hat{\rho}_h^{n}\mathcal U(\hat{\rho}_h^{n},\hat{s}_h^{n+1})-\hat{\rho}_h^{n}\mathcal U(\hat{\rho}_h^{n},\hat{s}_h^{n})}{\hat{s}_h^{n+1}-\hat{s}_h^n} \right) = 0 ~ , \\[2mm] &\frac{\mathbf s^{n+1}- \mathbf s^n}{\Delta t} + \mathbb D \hat{\Pi}^{2}[\hat{s_h^{n}} \vec{\boldsymbol \Lambda}^v] \mathbf u^{n+1/2} = 0 ~ , \\[2mm] \end{align}\end{split}\]

where \(\hat{l}^3(f)\) denotes the vector representing the linear form \(v_h \mapsto \int_{\hat{\Omega}} f(\boldsymbol \eta) v_h(\boldsymbol \eta) d \boldsymbol \eta\), that is the vector with components

\[\hat{l}^3(f)_{ijk}=\int_{\hat{\Omega}} f \Lambda^3_{ijk} \textrm d \boldsymbol \eta\]

and the weights in the the BasisProjectionOperator and the WeightedMassOperator are given by

\[\hat{\mathbf{u}}_h^{k} = (\mathbf{u}^{k})^\top \vec{\boldsymbol \Lambda}^v \in (V_h^0)^3 \, \text{for k in} \{n, n+1/2, n+1\}, \qquad \hat{s}_h^{k} = (s^{k})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \, \text{for k in} \{n, n+1/2, n+1\} \qquad \hat{\rho}_h^{n} = (\rho^{n})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \.\]

class Variables

Bases: object

property s: FEECVariable

property u: FEECVariable

class Options(model: Literal['full'] = 'full', gamma: float = 1.6666666666666667, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None, nonlin_solver: struphy.linear_algebra.solver.NonlinearSolverParameters = None, rho: struphy.models.variables.FEECVariable = None)

Bases: object

OptsModel

alias of Literal[‘full’]

model: Literal['full'] = 'full'

gamma: float = 1.6666666666666667

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

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

solver_params: SolverParameters = None

nonlin_solver: NonlinearSolverParameters = None

rho: FEECVariable = None

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_fields.VariationalMagFieldEvolve

Bases: Propagator

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

\[\begin{split}&\int_\Omega \partial_t (\rho \mathbf u) \cdot \mathbf v\,\textrm d \mathbf x - \int_\Omega \mathbf B \cdot \nabla \times (\mathbf \tilde{B} \times \mathbf v) \,\textrm d \mathbf x = 0 \qquad \forall \, \mathbf v \in (H^1)^3\,, \\[4mm] &\partial_t \mathbf B + \nabla \cdot ( \mathbf \tilde{B} \times \mathbf u ) = 0 \,.\end{split}\]

Where \(\tilde{\mathbf B}\) is either \(\mathbf B\) for full-f models, \(\mathbf B_0\) for linear models or \(\mathbf B_0+\mathbf B\) for \(\delta f\) models.

On the logical domain:

\[\begin{split}\begin{align} &\int_{\hat{\Omega}} \partial_t ( \hat{\rho}^3 \hat{\mathbf{u}}) \cdot G \hat{\mathbf{v}} \, \textrm d \boldsymbol \eta - \int_{\hat{\Omega}} \hat{\mathbf{B}}^2 \cdot G \,\nabla \times (\hat{\mathbf{B}}^2 \times \hat{\mathbf{v}}) \,\frac{1}{\sqrt g}\, \textrm d \boldsymbol \eta = 0 ~ , \\[2mm] &\partial_t \hat{\mathbf{B}}^2 + \nabla \times (\hat{\mathbf{B}}^2 \times \hat{\mathbf{u}}) = 0 ~ . \end{align}\end{split}\]

It is discretized as

\[\begin{split}\begin{align} &\mathbb M^v[\hat{\rho}_h^{n}] \frac{ \mathbf u^{n+1}-\mathbf u^n}{\Delta t} - (\mathbb C \hat{\Pi}^{1}[B_h^{n+1}} \cdot \vec{\boldsymbol \Lambda}^v])^\top \mathbb M^2 B^{n+\frac{1}{2}} \big) = 0 ~ , \\[2mm] &\frac{\mathbf b^{n+1}- \mathbf b^n}{\Delta t} + \mathbb C \hat{\Pi}^{1}[\hat{B_h^{n}} \cdot \vec{\boldsymbol \Lambda}^v]] \mathbf u^{n+1/2} = 0 ~ , \end{align}\end{split}\]

where weights in the the BasisProjectionOperator and the WeightedMassOperator are given by

\[\hat{\mathbf{B}}_h^{n+1/2} = (\mathbf{b}^{n+1/2})^\top \vec{\boldsymbol \Lambda}^2 \in V_h^2 \, \qquad \hat{\rho}_h^{n} = (\boldsymbol \rho^{n})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \,.\]

class Variables

Bases: object

property u: FEECVariable

property b: FEECVariable

class Options(model: Literal['full', 'full_p', 'linear'] = 'full', solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None, nonlin_solver: struphy.linear_algebra.solver.NonlinearSolverParameters = None)

Bases: object

OptsModel

alias of Literal[‘full’, ‘full_p’, ‘linear’]

model: Literal['full', 'full_p', 'linear'] = 'full'

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

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

solver_params: SolverParameters = None

nonlin_solver: NonlinearSolverParameters = None

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_fields.VariationalPBEvolve

Bases: Propagator

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

\[\begin{split}&\int_\Omega \partial_t (\rho \mathbf u) \cdot \mathbf v\,\textrm d \mathbf x - \int_\Omega \mathbf B \cdot \nabla \times (\tilde{\mathbf B} \times \mathbf v) - \int_\Omega \frac{1}{\gamma -1} (\nabla \cdot (\tilde{p} \mathbf v))\,\textrm d \mathbf x = 0 \qquad \forall \, \mathbf v \in (H^1)^3\,, \\[4mm] &\partial_t \mathbf B + \nabla \cdot ( \tilde{\mathbf B} \times \mathbf u ) = 0 \,, \\[4mm] &\partial_t p + \nabla \cdot(\tilde{p} \mathbf u) + (\gamma - 1) \tilde{p} \nabla \cdot u = 0 \,.\end{split}\]

Where \(\tilde{\mathbf B}\) (resp. \(\tilde{p}\)) is either \(\mathbf B\) (resp. \(p\)) for full-f models, \(\mathbf B_0\) (resp. \(p_0\)) for linear models or \(\mathbf B_0+\mathbf B\) (resp. \(p_0+p\)) for \(\delta f\) models.

On the logical domain:

\[\begin{split}\begin{align} &\int_{\hat{\Omega}} \partial_t ( \hat{\rho}^3 \hat{\mathbf{u}}) \cdot G \hat{\mathbf{v}} \, \textrm d \boldsymbol \eta - \int_{\hat{\Omega}} \hat{\mathbf{B}}^2 \cdot G \,\nabla \times (\hat{\mathbf{B}}^2 \times \hat{\mathbf{v}}) \,\frac{1}{\sqrt g}\, - \frac{g}{\gamma -1} \nabla \cdot (\hat{p} \hat{v}) \textrm d \boldsymbol \eta = 0 ~ , \\[2mm] &\partial_t \hat{\mathbf{B}}^2 + \nabla \times (\hat{\mathbf{B}}^2 \times \hat{\mathbf{u}}) = 0 ~ , \\[2mm] &\partial_t \hat{p} + \nabla \cdot (\hat{p} \hat{u}) + (\gamma - 1 ) \hat{p} \nabla \cdot G \hat{u} = 0 \, \end{align}\end{split}\]

It is discretized as

\[\begin{split}\begin{align} &\mathbb M^v[\hat{\rho}_h^{n}] \frac{ \mathbf u^{n+1}-\mathbf u^n}{\Delta t} - (\mathbb C \hat{\Pi}^{1}[\hat{\mathbf B_h^{n+\frac{1}{2}}} \cdot \vec{\boldsymbol \Lambda}^v])^\top \mathbb M^2 \mathbf B^{n+\frac{1}{2}} = 0 ~ , - (\mathbb D \hat{\Pi}^{2}[\hat{p_h^{n+\frac{1}{2}}} \cdot \vec{\boldsymbol \Lambda}^v])^\top \hat{l}^3(\frac{g}{\gamma-1}) \\[2mm] &\frac{\mathbf b^{n+1}- \mathbf b^n}{\Delta t} + \mathbb C \hat{\Pi}^{1}[\hat{\mathbf B_h^{n+\frac{1}{2}}} \cdot \vec{\boldsymbol \Lambda}^v]] \mathbf u^{n+1/2} = 0 ~ , \\[2mm] &\frac{\mathbf p^{n+1}- \mathbf p^n}{\Delta t} + \big(\mathbb D \hat{\Pi}^{2}[\hat{p_h^{n+\frac{1}{2}}} \cdot \vec{\boldsymbol \Lambda}^v]] + (\gamma - 1)\hat{\Pi}^{3}[\hat{p_h^{n+\frac{1}{2}}} \cdot \vec{\boldsymbol \Lambda}^3] \mathbb D \mathcal{U}^v \big) \mathbf u^{n+1/2}= 0 ~ , \\[2mm] \end{align}\end{split}\]

with

\[\hat{l}^3(f)_{ijk}=\int_{\hat{\Omega}} f \Lambda^3_{ijk} \textrm d \boldsymbol \eta\]

where weights in the the BasisProjectionOperator and the WeightedMassOperator are given by

\[\hat{\mathbf{B}}_h^{n+1/2} = (\mathbf{b}^{n+\frac{1}{2}})^\top \vec{\boldsymbol \Lambda}^2 \in V_h^2 \, \qquad \hat{\rho}_h^{n} = (\boldsymbol \rho^{n})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \, \qquad \hat{p}_h^{n+1/2} = (\boldsymbol p^{n+1/2})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \,.\]

and \(\mathcal{U}^v\) is BasisProjectionOperators.

class Variables

Bases: object

property p: FEECVariable

property u: FEECVariable

property b: FEECVariable

class Options(model: Literal['full_p', 'linear', 'deltaf'] = 'full_p', gamma: float = 1.6666666666666667, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None, nonlin_solver: struphy.linear_algebra.solver.NonlinearSolverParameters = None, div_u: struphy.models.variables.FEECVariable = None, u2: struphy.models.variables.FEECVariable = None, pt3: struphy.models.variables.FEECVariable = None, bt2: struphy.models.variables.FEECVariable = None)

Bases: object

OptsModel

alias of Literal[‘full_p’, ‘linear’, ‘deltaf’]

model: Literal['full_p', 'linear', 'deltaf'] = 'full_p'

gamma: float = 1.6666666666666667

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

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

solver_params: SolverParameters = None

nonlin_solver: NonlinearSolverParameters = None

div_u: FEECVariable = None

u2: FEECVariable = None

pt3: FEECVariable = None

bt2: FEECVariable = None

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_fields.VariationalQBEvolve

Bases: Propagator

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

\[\begin{split}&\int_\Omega \partial_t (\rho \mathbf u) \cdot \mathbf v\,\textrm d \mathbf x - \int_\Omega \mathbf B \cdot \nabla \times (\tilde{\mathbf B} \times \mathbf v) - \int_\Omega \frac{2 q}{\gamma -1} (\nabla \cdot (\tilde{q} \mathbf v))\,\textrm d \mathbf x = 0 \qquad \forall \, \mathbf v \in (H^1)^3\,, \\[4mm] &\partial_t \mathbf B + \nabla \cdot ( \tilde{\mathbf B} \times \mathbf u ) = 0 \,, \\[4mm] &\partial_t q + \nabla \cdot(\tilde{q} \mathbf u) + (\gamma/2 - 1) \tilde{q} \nabla \cdot u = 0 \,.\end{split}\]

Where \(\tilde{\mathbf B}\) (resp. \(\tilde{q}\)) is either \(\mathbf B\) (resp. \(q\)) for full-f models, \(\mathbf B_0\) (resp. \(q_0\)) for linear models or \(\mathbf B_0+\mathbf B\) (resp. \(q_0+q\)) for \(\delta f\) models.

On the logical domain:

\[\begin{split}\begin{align} &\int_{\hat{\Omega}} \partial_t ( \hat{\rho}^3 \hat{\mathbf{u}}) \cdot G \hat{\mathbf{v}} \, \textrm d \boldsymbol \eta - \int_{\hat{\Omega}} \hat{\mathbf{B}}^2 \cdot G \,\nabla \times (\hat{\mathbf{B}}^2 \times \hat{\mathbf{v}}) \,\frac{1}{\sqrt g}\, - \frac{q}{\gamma -1} \nabla \cdot (\hat{q} \hat{v}) \textrm d \boldsymbol \eta = 0 ~ , \\[2mm] &\partial_t \hat{\mathbf{B}}^2 + \nabla \times (\hat{\mathbf{B}}^2 \times \hat{\mathbf{u}}) = 0 ~ , \\[2mm] &\partial_t \hat{q} + \nabla \cdot (\hat{q} \hat{u}) + (\gamma/2 - 1 ) \hat{p} \nabla \cdot G \hat{u} = 0 \, \end{align}\end{split}\]

It is discretized as

\[\begin{split}\begin{align} &\mathbb M^v[\hat{\rho}_h^{n}] \frac{ \mathbf u^{n+1}-\mathbf u^n}{\Delta t} - (\mathbb C \hat{\Pi}^{1}[\hat{\mathbf B_h^{n}} \cdot \vec{\boldsymbol \Lambda}^v])^\top \mathbb M^2 \mathbf B^{n+\frac{1}{2}}- \Big(\big(\mathbb D \hat{\Pi}^{2}[\hat{q_h^n} \cdot \vec{\boldsymbol \Lambda}^v]] + (\gamma/2 - 1)\hat{\Pi}^{3}[\hat{q_h^n} \cdot \vec{\boldsymbol \Lambda}^3] \mathbb D \mathcal{U}^v \big) \mathbf v \Big)^\top \hat{l}^3(\frac{2q^{n+\frac{1}{2}}}{\gamma-1}) = 0 ~ , \\[2mm] &\frac{\mathbf b^{n+1}- \mathbf b^n}{\Delta t} + \mathbb C \hat{\Pi}^{1}[\hat{\mathbf B_h^{n+\frac{1}{2}}} \cdot \vec{\boldsymbol \Lambda}^v]] \mathbf u^{n+1/2} = 0 ~ , \\[2mm] &\frac{\mathbf q^{n+1}- \mathbf q^n}{\Delta t} + \big(\mathbb D \hat{\Pi}^{2}[\hat{q_h^n} \cdot \vec{\boldsymbol \Lambda}^v]] + (\gamma/2 - 1)\hat{\Pi}^{3}[\hat{q_h^n} \cdot \vec{\boldsymbol \Lambda}^3] \mathbb D \mathcal{U}^v \big) \mathbf u^{n+1/2}= 0 ~ , \\[2mm] \end{align}\end{split}\]

with

\[\hat{l}^3(f)_{ijk}=\int_{\hat{\Omega}} f \Lambda^3_{ijk} \textrm d \boldsymbol \eta\]

where weights in the the BasisProjectionOperator and the WeightedMassOperator are given by

\[\hat{\mathbf{B}}_h^{n+1/2} = (\mathbf{b}^{n+\frac{1}{2}})^\top \vec{\boldsymbol \Lambda}^2 \in V_h^2 \, \qquad \hat{\rho}_h^{n} = (\boldsymbol \rho^{n})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \, \qquad \hat{q}_h^{n+1/2} = (\boldsymbol q^{n+1/2})^\top \vec{\boldsymbol \Lambda}^3 \in V_h^3 \,.\]

and \(\mathcal{U}^v\) is BasisProjectionOperators.

class Variables

Bases: object

property q: FEECVariable

property u: FEECVariable

property b: FEECVariable

class Options(model: Literal['full_q', 'linear_q', 'deltaf_q'] = 'full_q', gamma: float = 1.6666666666666667, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None, nonlin_solver: struphy.linear_algebra.solver.NonlinearSolverParameters = None, div_u: struphy.models.variables.FEECVariable = None, u2: struphy.models.variables.FEECVariable = None, qt3: struphy.models.variables.FEECVariable = None, bt2: struphy.models.variables.FEECVariable = None)

Bases: object

OptsModel

alias of Literal[‘full_q’, ‘linear_q’, ‘deltaf_q’]

model: Literal['full_q', 'linear_q', 'deltaf_q'] = 'full_q'

gamma: float = 1.6666666666666667

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

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

solver_params: SolverParameters = None

nonlin_solver: NonlinearSolverParameters = None

div_u: FEECVariable = None

u2: FEECVariable = None

qt3: FEECVariable = None

bt2: FEECVariable = None

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_fields.VariationalViscosity

Bases: Propagator

FEEC discretization of the following equations: find \(s \in L^2\) and \(\mathbf u \in (H^1)^3\) such that

\[\begin{split}&\int_\Omega \partial_t (\rho \mathbf u) \cdot \mathbf v\,\textrm d \mathbf x + \int_\Omega (\mu + \mu_a(\mathbf x)) \nabla \mathbf u : \nabla \mathbf v \,\textrm d \mathbf x = 0 \qquad \forall \, \mathbf v \in (H^1)^3 \,, \\[4mm] &\int_\Omega \frac{\partial \mathcal U}{\partial s} \partial_t s \, q \,\textrm d \mathbf x - \mu \int_\Omega |\nabla \mathbf u|^2 \, q \,\textrm d \mathbf x = 0 \qquad \forall \, q \in L^2\,\text{if using } s, \\[4mm] &\int_\Omega \frac{1}{\gamma - 1} \partial_t p \, q\,\textrm d \mathbf x - \mu \int_\Omega |\nabla \mathbf u|^2 \, q \,\textrm d \mathbf x = 0 \qquad \forall \, q \in L^2\, \text{if using } p.\end{split}\]

With \(\mu_a(\mathbf x) = \mu_a |\nabla \mathbf u(\mathbf x)|\)

On the logical domain:

\[\begin{split}\begin{align} &\int_{\hat{\Omega}} \partial_t ( \hat{\rho}^3 \hat{\mathbf{u}}) \cdot G \hat{\mathbf{v}} \, \textrm d \boldsymbol \eta + \mu \int_{\hat{\Omega}} \nabla (DF \hat{\mathbf{u}}) : \nabla (DF \hat{\mathbf{v}}) \,\frac{1}{\sqrt g}\, \textrm d \boldsymbol \eta = 0 ~ , \\[2mm] &\int_{\hat{\Omega}} \partial_t (\hat{\rho} \hat{e}(\hat{\rho}, \hat{s})) \hat{w} \,\frac{1}{\sqrt g}\, \textrm d \boldsymbol \eta - \int_{\hat{\Omega}} (\mu + \mu_a(\boldsymbol \eta)) \nabla (DF \hat{\mathbf{u}}) : \nabla (DF \hat{\mathbf{u}}) \hat{w} \, \textrm d \boldsymbol \eta = 0 ~ , \text{if using } s, \\[2mm] &\int_{\hat{\Omega}} \partial_t (\frac{1}{\gamma -1} \hat{p} ) \hat{w} \,\frac{1}{\sqrt g}\, \textrm d \boldsymbol \eta - \int_{\hat{\Omega}} (\mu + \mu_a(\boldsymbol \eta)) \nabla (DF \hat{\mathbf{u}}) : \nabla (DF \hat{\mathbf{u}}) \hat{w} \, \textrm d \boldsymbol \eta = 0 ~, \text{if using } p. \end{align}\end{split}\]

It is discretized as

\[\begin{split}\begin{align} &\mathbb M^v[\hat{\rho}_h^{n}] \frac{ \mathbf u^{n+1}-\mathbf u^n}{\Delta t} + \sum_\nu (\mathbb G \mathcal{X}^v_\nu)^T (\mu \mathbb M_0 + \mu_a \mathbb M_0[|\nabla u|] \mathbb G \mathcal{X}^v_\nu \mathbf u^{n+1} = 0 ~ , \\[2mm] &\frac{P^{3}(\hat{\rho}_h^{n}\mathcal U(\hat{\rho}_h^{n},\hat{s}_h^{n}))- P^{3}(\hat{\rho}_h^{n}\mathcal U(\hat{\rho}_h^{n},\hat{s}_h^{n+1}))}{\Delta t} - \mu P^3(\sum_\nu DF \mathcal{X}^v_\nu \frac{ \mathbf u^{n+1}+\mathbf u^n}{2} \cdot DF \mathcal{X}^v_\nu \mathbf u^{n+1}) = 0 ~ , \text{if using } s, \\[2mm] &\frac{1}{\gamma -1}\frac{p^{n+1}- p^{n}}{\Delta t} - \mu P^3(\sum_\nu DF \mathcal{X}^v_\nu \frac{ \mathbf u^{n+1}+\mathbf u^n}{2} \cdot DF \mathcal{X}^v_\nu \mathbf u^{n+1}) = 0 ~ , \text{if using } p. \end{align}\end{split}\]

where $P^3$ denotes the $L^2$ projection in the last space of the de Rham sequence and the weights in \(\mathbb M_0[|\nabla u|]\) are given by

\[P^0(g \sqrt{\sum_\nu |(\mathbb G \mathcal{X}^v_\nu \mathbb u)^\top \vec{\boldsymbol \Lambda}^0 |^2]})^\top \vec{\boldsymbol \Lambda}^0 ~.\]

class Variables

Bases: object

property s: FEECVariable

property u: FEECVariable

class Options(model: Literal['full', 'full_p', 'full_q', 'linear_p', 'linear_q', 'deltaf_q'] = 'full', gamma: float = 1.6666666666666667, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixDiagonalPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None, nonlin_solver: struphy.linear_algebra.solver.NonlinearSolverParameters = None, rho: struphy.models.variables.FEECVariable = None, pt3: struphy.models.variables.FEECVariable = None, mu: float = 0.0, mu_a: float = 0.0, alpha: float = 0.0)

Bases: object

OptsModel

alias of Literal[‘full’, ‘full_p’, ‘full_q’, ‘linear_p’, ‘linear_q’, ‘deltaf_q’]

model: Literal['full', 'full_p', 'full_q', 'linear_p', 'linear_q', 'deltaf_q'] = 'full'

gamma: float = 1.6666666666666667

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

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

solver_params: SolverParameters = None

nonlin_solver: NonlinearSolverParameters = None

rho: FEECVariable = None

pt3: FEECVariable = None

mu: float = 0.0

mu_a: float = 0.0

alpha: float = 0.0

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_fields.VariationalResistivity

Bases: Propagator

FEEC discretization of the following equations: find \(s \in L^2\) and \(\mathbf B \in H(\textrm{div})\) such that

\[\begin{split}&\int_\Omega \partial_t \mathbf B \cdot \mathbf v \,\textrm d \mathbf x + \int_\Omega (\eta + \eta_a(\mathbf x)) \nabla \times \mathbf B \cdot \nabla \times \mathbf v \,\textrm d \mathbf x = 0 \qquad \forall \, \mathbf v \in H(\textrm{div}) \,, \\[4mm] &\int_\Omega \frac{\partial \mathcal U}{\partial s} \partial_t s \, q\,\textrm d \mathbf x - \int_\Omega (\eta + \eta_a(\mathbf x)) |\nabla \times \mathbf B|^2 \, q \,\textrm d \mathbf x = 0 \qquad \forall \, q \in L^2\, \text{if using } s, \\[4mm] &\int_\Omega \frac{1}{\gamma - 1} \partial_t p \, q\,\textrm d \mathbf x - \int_\Omega (\eta + \eta_a(\mathbf x)) |\nabla \times \mathbf B|^2 \, q \,\textrm d \mathbf x = 0 \qquad \forall \, q \in L^2\, \text{if using } p.\end{split}\]

With \(\eta_a(\mathbf x) = \eta_a |\nabla \times \mathbf B(\mathbf x)|\)

On the logical domain:

\[\begin{split}\begin{align} &\partial_t \hat{\boldsymbol B} - \eta \Delta \hat{\boldsymbol B} = 0 ~ , \\[2mm] &\int_{\hat{\Omega}} \partial_t (\hat{\rho} \mathcal U(\hat{\rho}, \hat{s})) \hat{w} \,\frac{1}{\sqrt g}\, \textrm d \boldsymbol \eta - \int_{\hat{\Omega}} (\eta + \eta_a(\mathbf \eta)) |DF^{-T}\tilde{\nabla} \times \hat{\boldsymbol B}|^2 \hat{w} \, \textrm d \boldsymbol \eta = 0 ~ , \text{if using } s, \\[2mm] &\int_{\hat{\Omega}} \partial_t (\frac{1}{\gamma -1} \hat{p} ) \hat{w} \,\frac{1}{\sqrt g}\, \textrm d \boldsymbol \eta - \int_{\hat{\Omega}} (\eta + \eta_a(\mathbf \eta)) |DF^{-T}\tilde{\nabla} \times \hat{\boldsymbol B}|^2 \hat{w} \, \textrm d \boldsymbol \eta = 0 ~, \text{if using } p. \end{align}\end{split}\]

It is discretized as

\[\begin{split}\begin{align} &\frac{\mathbf B^{n+1}-\mathbf B^n}{\Delta t} + \, \mathbb C \mathbb M_1^{-1} (\eta M_1 + \eta_a M_1[|\nabla \times \mathbf B|]) M_1^{-1} \mathbb C^T \mathbb M_2 \mathbf B^{n+1} = 0 ~ , \\[2mm] &\frac{P^{3}(\rho e(s^{n+1})- P^{3}(\rho e(s^{n}))}{\Delta t} - P^3((\eta + \eta_a(\mathbf x)) DF^{-T} \tilde{\mathbb C} \frac{ \mathbf B^{n+1}+\mathbf B^n}{2} \cdot DF^{-T} \tilde{\mathbb C} \mathbf B^{n+1}) = 0 ~ , \text{if using } s, \\[2mm] &\frac{1}{\gamma -1}\frac{p^{n+1}- p^{n}}{\Delta t} - P^3((\eta + \eta_a(\mathbf x)) DF^{-T} \tilde{\mathbb C} \frac{ \mathbf B^{n+1}+\mathbf B^n}{2} \cdot DF^{-T} \tilde{\mathbb C} \mathbf B^{n+1}) = 0 ~ , \text{if using } p. \end{align}\end{split}\]

where $P^3$ denotes the $L^2$ projection in the last space of the de Rham sequence.

class Variables

Bases: object

property s: FEECVariable

property b: FEECVariable

class Options(model: Literal['full', 'full_p', 'full_q', 'linear_p', 'linear_q', 'deltaf_q'] = 'full', gamma: float = 1.6666666666666667, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixDiagonalPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None, nonlin_solver: struphy.linear_algebra.solver.NonlinearSolverParameters = None, linearize_current: bool = False, rho: struphy.models.variables.FEECVariable = None, pt3: struphy.models.variables.FEECVariable = None, eta: float = 0.0, eta_a: float = 0.0)

Bases: object

OptsModel

alias of Literal[‘full’, ‘full_p’, ‘full_q’, ‘linear_p’, ‘linear_q’, ‘deltaf_q’]

model: Literal['full', 'full_p', 'full_q', 'linear_p', 'linear_q', 'deltaf_q'] = 'full'

gamma: float = 1.6666666666666667

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

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

solver_params: SolverParameters = None

nonlin_solver: NonlinearSolverParameters = None

linearize_current: bool = False

rho: FEECVariable = None

pt3: FEECVariable = None

eta: float = 0.0

eta_a: float = 0.0

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_fields.TimeDependentSource

Bases: Propagator

Propagates a source term \(S(t) \in V_h^n\) of the form

\[S(t) = \sum_{ijk} c_{ijk} \Lambda^n_{ijk} * h(\omega t)\,,\]

where \(h(\omega t)\) is one of the functions in Notes.

Notes

• \(h(\omega t) = \cos(\omega t)\) (default)

• \(h(\omega t) = \sin(\omega t)\)

class Variables

Bases: object

property source: FEECVariable

class Options(omega: float = 6.283185307179586, hfun: Literal['cos', 'sin'] = 'cos')

Bases: object

OptsTimeSource

alias of Literal[‘cos’, ‘sin’]

omega: float = 6.283185307179586

hfun: Literal['cos', 'sin'] = 'cos'

allocate()

Allocate all data/objects of the instance.

class struphy.propagators.propagators_fields.AdiabaticPhi(phi: StencilVector, *, A_mat: WeightedMassOperator = 'M0', rho: StencilVector | tuple | None = None, sigma_1: float = 1.0, x0: StencilVector | None = None, **params)

Bases: Propagator

Electrostatic potential for adiabatic electrons, computed from

\[n_e = n_{e0}\,\exp \left( \frac{e \phi}{k_B T_e} \right) \approx n_{e0} \left( 1 + \frac{e \phi}{k_B T_{e0}} \right)\,,\]

where \(n_{e0}\) and \(T_{e0}\) denote electron equilibrium density and temperature, respectively.

This is solved in weak form: find \(\phi \in H^1\) such that

\[\int_\Omega \psi\, \frac{n_{e0}(\mathbf x)}{T_{e0}(\mathbf x)}\phi \,\textrm d \mathbf x = \int_\Omega \psi\, (n_{e}(\mathbf x) - n_{e0}(\mathbf x))\,\textrm d \mathbf x \qquad \forall \ \psi \in H^1\,.\]

The equation is discretized as

\[\sigma_1 \mathbb M^0_{n/T} \boldsymbol \phi = (\Lambda^0, n_{e} - n_{e0} )_{L^2}\,,\]

where \(M^0_{n/T}\) is a WeightedMassOperator and \(\sigma_1\) is a normalization parameter.

Parameters:

• phi (StencilVector) – FE coefficients of the solution as a discrete 0-form.

• A_mat (WeightedMassOperator) – The matrix to invert.

• rho (StencilVector or tuple) – 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. In case b) the first tuple entry must be AccumulatorVector, and the second entry must be Particles.

• sigma_1 (float) – Normalization parameter.

• x0 (StencilVector) – Initial guess for the iterative solver (optional, can be set with a setter later).

• **params (dict) – Parameters for the iterative solver (see __init__ for details).

property rho

Right-hand side FE coefficients of a 0-form. Can be either a) StencilVector or b) 2-tuple. In the latter case, the first tuple entry must be AccumulatorVector, and the second entry must be Particles.

property x0

psydac.linalg.stencil.StencilVector or struphy.polar.basic.PolarVector. First guess of the iterative solver.

class struphy.propagators.propagators_fields.HasegawaWakatani

Bases: Propagator

FEEC discretization of the following equations: find \((n, \omega) \in H^1 \times H^1\) such that

\[\begin{split}&\int_\Omega\frac{\partial n}{\partial t} m \,\textrm d \mathbf x = \int_\Omega C(x, y)(\phi - n) \, m \,\textrm d \mathbf x - \int_\Omega \phi [n, m] \,\textrm d \mathbf x - \kappa \int_\Omega \partial_y \phi \,m \,\textrm d \mathbf x - \nu \int_\Omega \nabla n \cdot \nabla m \,\textrm d \mathbf x \qquad \forall m \in H^1\,, \\[2mm] &\int_\Omega\frac{\partial \omega}{\partial t} \psi \,\textrm d \mathbf x = \int_\Omega C(x, y)(\phi - n) \, \psi \,\textrm d \mathbf x - \int_\Omega \phi [\omega, \psi] \,\textrm d \mathbf x - \nu \int_\Omega \nabla \omega \cdot \nabla \psi \,\textrm d \mathbf x \qquad \forall \psi \in H^1\,,\end{split}\]

where \(\phi \in H^1\) is a given stream function, \(C = C(x, y)\), \(\kappa\) and \(\nu\) are constants and \([a, b] = \partial_x a \partial_y b - \partial_y a \partial_x b\).

Time discretization: explicit Runge-Kutta, see ODEsolverFEEC.

Parameters:

• n0 (StencilVector) – The density.

• omega0 (StencilVector) – The stream function.

• phi (SplineFuncion) – The potential.

• c_fun (str) – Defines the function c(x,y) in front of (phi - n).

• kappa (float) – Equation parameters.

• nu (float) – Equation parameters.

• algo (str) – See ButcherTableau for available algorithms.

• M0_solver (dict) – Solver parameters for M0 inversion.

class Variables

Bases: object

property n: FEECVariable

property omega: FEECVariable

class Options(phi: struphy.models.variables.FEECVariable = None, c_fun: Literal['const'] = 'const', kappa: float = 1.0, nu: float = 0.01, butcher: struphy.ode.utils.ButcherTableau = None, solver: Literal['pcg', 'cg'] = 'pcg', precond: Literal['MassMatrixPreconditioner', 'MassMatrixDiagonalPreconditioner', None] = 'MassMatrixPreconditioner', solver_params: struphy.linear_algebra.solver.SolverParameters = None)

Bases: object

OptsCfun

alias of Literal[‘const’]

phi: FEECVariable = None

c_fun: Literal['const'] = 'const'

kappa: float = 1.0

nu: float = 0.01

butcher: ButcherTableau = None

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_fields.TwoFluidQuasiNeutralFull

Bases: Propagator

FEEC discretization of the following equations: find \(\mathbf u \in H(\textnormal{div})\), \(\mathbf u_e \in H(\textnormal{div})\) and \(\mathbf \phi \in L^2\) such that

\[\begin{split}\int_{\Omega} \partial_t \mathbf{u}\cdot \mathbf{v} \, \textrm d\mathbf{x} &= \int_{\Omega} \phi \nabla \! \cdot \! \mathbf{v} \, \textrm d\mathbf{x} + \int_{\Omega} \mathbf{u}\! \times \! \mathbf{B}_0 \cdot \mathbf{v} \, \textrm d\mathbf{x} + \nu \int_{\Omega} \nabla \mathbf{u}\! : \! \nabla \mathbf{v} \, \textrm d\mathbf{x} + \int_{\Omega} f \mathbf{v} \, \textrm d\mathbf{x} \qquad \forall \, \mathbf{v} \in H(\textrm{div}) \,. \\[2mm] 0 &= - \int_{\Omega} \phi \nabla \! \cdot \! \mathbf{v_e} \, \textrm d\mathbf{x} - \int_{\Omega} \mathbf{u_e} \! \times \! \mathbf{B}_0 \cdot \mathbf{v_e} \, \textrm d\mathbf{x} + \nu_e \int_{\Omega} \nabla \mathbf{u_e} \!: \! \nabla \mathbf{v_e} \, \textrm d\mathbf{x} + \int_{\Omega} f_e \mathbf{v_e} \, \textrm d\mathbf{x} \qquad \forall \ \mathbf{v_e} \in H(\textrm{div}) \,. \\[2mm] 0 &= \int_{\Omega} \psi \nabla \cdot (\mathbf{u}-\mathbf{u_e}) \, \textrm d\mathbf{x} \qquad \forall \, \psi \in L^2 \,.\end{split}\]

Time discretization: fully implicit.

class Variables

Bases: object

property u: FEECVariable

property ue: FEECVariable

property phi: FEECVariable

class Options(nu: float = 1.0, nu_e: float = 0.01, eps_norm: float = 1.0, solver: Literal['pbicgstab', 'bicgstab', 'GMRES'] = 'GMRES', solver_params: struphy.linear_algebra.solver.SolverParameters = None, a: float = 1.0, R0: float = 1.0, B0: float = 10.0, Bp: float = 12.0, alpha: float = 0.1, beta: float = 1.0, stab_sigma: float = 1e-05, variant: Literal['Uzawa', 'GMRES'] = 'Uzawa', method_to_solve: Literal['SparseSolver', 'ScipySparse', 'InexactNPInverse', 'DirectNPInverse'] = 'DirectNPInverse', preconditioner: bool = False, spectralanalysis: bool = False, lifting: bool = False, dimension: Literal['1D', '2D', 'Restelli', 'Tokamak'] = '2D', D1_dt: float = 0.001)

Bases: object

OptsDimension

alias of Literal[‘1D’, ‘2D’, ‘Restelli’, ‘Tokamak’]

nu: float = 1.0

nu_e: float = 0.01

eps_norm: float = 1.0

solver: Literal['pbicgstab', 'bicgstab', 'GMRES'] = 'GMRES'

solver_params: SolverParameters = None

a: float = 1.0

R0: float = 1.0

B0: float = 10.0

Bp: float = 12.0

alpha: float = 0.1

beta: float = 1.0

stab_sigma: float = 1e-05

variant: Literal['Uzawa', 'GMRES'] = 'Uzawa'

method_to_solve: Literal['SparseSolver', 'ScipySparse', 'InexactNPInverse', 'DirectNPInverse'] = 'DirectNPInverse'

preconditioner: bool = False

spectralanalysis: bool = False

lifting: bool = False

dimension: Literal['1D', '2D', 'Restelli', 'Tokamak'] = '2D'

D1_dt: float = 0.001

allocate()

Allocate all data/objects of the instance.