Weighted mass operators#

class struphy.feec.mass.WeightedMassOperators(derham: Derham, domain: Domain, matrix_free: bool = False, verbose: bool = True, **weights)[source]#

Bases: object

Collection of pre-defined struphy.feec.mass.WeightedMassOperator.

Parameters:

derham (Derham) – Discrete de Rham sequence on the logical unit cube.
domain (Geometry) – Mapping from logical unit cube to physical domain and corresponding metric coefficients.
**weights (dict) – Objects to access callables that can serve as weight functions.
matrix_free (bool) – If set to true will not compute the matrix associated with the operator but directly compute the product when called
verbose (bool) – Show info on screen.

Notes

Possible choices for key-value pairs in weights are, at the moment:

eq_mhd: MHDequilibrium

property derham#: Discrete de Rham sequence on the logical unit cube.

property domain#: Mapping from the logical unit cube to the physical domain with corresponding metric coefficients.

property weights#: Dictionary of objects that provide access to callables that can serve as weight functions.

property verbose#

show info on screen.

Type:: Bool

property selected_weight#: String identifying one key of “weigths”. This key is used when selecting weight functions.

property M0#: Mass matrix

\[\mathbb M^0_{ijk, mno} = \int \Lambda^0_{ijk}\, \Lambda^0_{mno} \sqrt g\, \textnormal d \boldsymbol\eta.\]

property M1#: Mass matrix

\[\mathbb M^1_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^1_{\mu,ijk}\, G^{-1}\, \vec{\Lambda}^1_{\nu, mno} \sqrt g\, \textnormal d \boldsymbol\eta.\]

property M2#: Mass matrix

\[\mathbb M^2_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^2_{\mu,ijk}\, G\, \vec{\Lambda}^2_{\nu, mno} \frac{1}{\sqrt g}\, \textnormal d \boldsymbol\eta.\]

property M3#: Mass matrix

\[\mathbb M^3_{ijk, mno} = \int \Lambda^3_{ijk}\, \Lambda^3_{mno} \frac{1}{\sqrt g}\, \textnormal d \boldsymbol\eta.\]

property Mv#: Mass matrix

\[\mathbb M^v_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^v_{\mu,ijk}\, G\, \vec{\Lambda}^v_{\nu, mno} \sqrt g\, \textnormal d \boldsymbol\eta.\]

property M1n#

Mass matrix

\[\mathbb M^{1,n}_{(\mu,ijk), (\nu,mno)} = \int n^0_{\textnormal{eq}}(\boldsymbol \eta) \vec{\Lambda}^1_{\mu,ijk}\, G^{-1}\, \vec{\Lambda}^1_{\nu, mno} \sqrt g\, \textnormal d \boldsymbol\eta.\]

where \(n^0_{\textnormal{eq}}(\boldsymbol \eta)\) is an MHD equilibrium density (0-form).

property M2n#

Mass matrix

\[\mathbb M^{2,n}_{(\mu,ijk), (\nu,mno)} = \int n^0_{\textnormal{eq}}(\boldsymbol \eta) \vec{\Lambda}^2_{\mu,ijk}\, G\, \vec{\Lambda}^2_{\nu, mno} \frac{1}{\sqrt g}\, \textnormal d \boldsymbol\eta.\]

where \(n^0_{\textnormal{eq}}(\boldsymbol \eta)\) is an MHD equilibrium density (0-form).

property Mvn#

Mass matrix

\[\mathbb M^{v,n}_{(\mu,ijk), (\nu,mno)} = \int n^0_{\textnormal{eq}}(\boldsymbol \eta) \vec{\Lambda}^v_{\mu,ijk}\, G\, \vec{\Lambda}^v_{\nu, mno} \sqrt g\, \textnormal d \boldsymbol\eta.\]

where \(n^0_{\textnormal{eq}}(\boldsymbol \eta)\) is an MHD equilibrium density (0-form).

property M1ninv#

Mass matrix

\[\mathbb M^{1,\frac{1}{n}}_{(\mu,ijk), (\nu,mno)} = \int \frac{1}{n^0_{\textnormal{eq}}(\boldsymbol \eta)} \vec{\Lambda}^1_{\mu,ijk}\, G^{-1}\, \vec{\Lambda}^1_{\nu, mno} \sqrt g\, \textnormal d \boldsymbol\eta.\]

where \(n^0_{\textnormal{eq}}(\boldsymbol \eta)\) is an MHD equilibrium density (0-form).

property M1J#

Mass matrix

\[\mathbb M^{1,J}_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^1_{\mu,ijk}\, G^{-1}\, \mathcal R^J\, \vec{\Lambda}^2_{\nu, mno} \, \textnormal d \boldsymbol\eta.\]

with the rotation matrix

\[\mathcal R^J_{\alpha, \nu} := \epsilon_{\alpha \beta \nu}\, J^2_{\textnormal{eq}, \beta}\,,\qquad s.t. \qquad \mathcal R^J \vec v = \vec J^2_{\textnormal{eq}} \times \vec v\,,\]

where \(\epsilon_{\alpha \beta \nu}\) stands for the Levi-Civita tensor and \(J^2_{\textnormal{eq}, \beta}\) is the \(\beta\)-component of the MHD equilibrium current density (2-form).

property M2J#

Mass matrix

\[\mathbb M^{2,J}_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^2_{\mu,ijk}\, \mathcal R^J\, \vec{\Lambda}^2_{\nu, mno} \, \frac{1}{\sqrt g}\, \textnormal d \boldsymbol\eta.\]

with the rotation matrix

\[\mathcal R^J_{\alpha, \nu} := \epsilon_{\alpha \beta \nu}\, J^2_{\textnormal{eq}, \beta}\,,\qquad s.t. \qquad \mathcal R^J \vec v = \vec J^2_{\textnormal{eq}} \times \vec v\,,\]

where \(\epsilon_{\alpha \beta \nu}\) stands for the Levi-Civita tensor and \(J^2_{\textnormal{eq}, \beta}\) is the \(\beta\)-component of the MHD equilibrium current density (2-form).

property MvJ#

Mass matrix

\[\mathbb M^{v,J}_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^v_{\mu,ijk}\, \mathcal R^J\, \vec{\Lambda}^v_{\nu, mno} \, \textnormal d \boldsymbol\eta.\]

with the rotation matrix

\[\mathcal R^J_{\alpha, \nu} := \epsilon_{\alpha \beta \nu}\, J^2_{\textnormal{eq}, \beta}\,,\qquad s.t. \qquad \mathcal R^J \vec v = \vec J^2_{\textnormal{eq}} \times \vec v\,,\]

where \(\epsilon_{\alpha \beta \nu}\) stands for the Levi-Civita tensor and \(J^2_{\textnormal{eq}, \beta}\) is the \(\beta\)-component of the MHD equilibrium current density (2-form).

property M2B_div0#

Mass matrix

\[\mathbb M^{2,B}_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^2_{\mu,ijk}\, \mathcal R^J\, \vec{\Lambda}^2_{\nu, mno} \, \frac{1}{\sqrt g}\, \textnormal d \boldsymbol\eta.\]

with the rotation matrix

\[\mathcal R^J_{\alpha, \nu} := \epsilon_{\alpha \beta \nu}\, B^2_{\textnormal{eq}, \beta}\,,\qquad s.t. \qquad \mathcal R^J \vec v = \vec B^2_{\textnormal{eq}} \times \vec v\,,\]

where \(\epsilon_{\alpha \beta \nu}\) stands for the Levi-Civita tensor and \(B^2_{\textnormal{eq}, \beta}\) is the \(\beta\)-component of the MHD equilibrium magnetic field (2-form).

property M2B#

Mass matrix

\[\mathbb M^{2,B}_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^2_{\mu,ijk}\, \mathcal R^J\, \vec{\Lambda}^2_{\nu, mno} \, \frac{1}{\sqrt g}\, \textnormal d \boldsymbol\eta.\]

with the rotation matrix

\[\mathcal R^J_{\alpha, \nu} := \epsilon_{\alpha \beta \nu}\, B^2_{\textnormal{eq}, \beta}\,,\qquad s.t. \qquad \mathcal R^J \vec v = \vec B^2_{\textnormal{eq}} \times \vec v\,,\]

where \(\epsilon_{\alpha \beta \nu}\) stands for the Levi-Civita tensor and \(B^2_{\textnormal{eq}, \beta}\) is the \(\beta\)-component of the MHD equilibrium magnetic field (2-form).

property M2Bn#

Mass matrix

\[\mathbb M^{2,BN}_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^2_{\mu,ijk}\, \mathcal R^J\, \vec{\Lambda}^2_{\nu, mno} \, \frac{1}{n^0_{\textnormal{eq}}(\boldsymbol \eta)}\, \frac{1}{\sqrt g}\, \textnormal d \boldsymbol\eta.\]

with the rotation matrix

\[\mathcal R^J_{\alpha, \nu} := \epsilon_{\alpha \beta \nu}\, B^2_{\textnormal{eq}, \beta}\,,\qquad s.t. \qquad \mathcal R^J \vec v = \vec B^2_{\textnormal{eq}} \times \vec v\,,\]

where \(\epsilon_{\alpha \beta \nu}\) stands for the Levi-Civita tensor and \(B^2_{\textnormal{eq}, \beta}\) is the \(\beta\)-component of the MHD equilibrium magnetic field (2-form).

property M1Bninv#

Mass matrix

\[\mathbb M^{1,B\frac{1}{n}}_{(\mu,ijk), (\nu,mno)} = \int \frac{1}{n^0_{\textnormal{eq}}(\boldsymbol \eta)}\, \vec{\Lambda}^1_{\mu,ijk}\, G^{-1}\, \mathcal R^J_{\alpha, \gamma}\, G^{-1}_{\gamma,\nu}\, \vec{\Lambda}^1_{\nu, mno} \, \sqrt g\, \textnormal d \boldsymbol\eta.\]

with the rotation matrix

\[\mathcal R^J_{\alpha, \nu} := \epsilon_{\alpha \beta \nu}\, B^2_{\textnormal{eq}, \beta}\,,\qquad s.t. \qquad \mathcal R^J \vec v = \vec B^2_{\textnormal{eq}} \times \vec v\,,\]

where \(\epsilon_{\alpha \beta \nu}\) stands for the Levi-Civita tensor and \(B^2_{\textnormal{eq}, \beta}\) is the \(\beta\)-component of the MHD equilibrium magnetic field (2-form).

property M1perp#: Mass matrix

\[\begin{split}\mathbb M^{1, \perp}_{(\mu, ijk), (\nu, mno)} = \int \vec{\Lambda}^1_{\mu, ijk}\, DF^{-1} \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix} DF^{-\top} \vec{\Lambda}^1_{\nu, mno} \sqrt g\, \textnormal d \boldsymbol\eta.\end{split}\]

property M0ad#: Mass matrix

\[\mathbb M^0_{ijk, mno} = \int \Lambda^0_{ijk}\, \Lambda^0_{mno} \sqrt g\, \textnormal d \boldsymbol\eta.\]

property M1gyro#

Mass matrix

\[\mathbb M^{1,n}_{(\mu,ijk), (\nu,mno)} = \int n^0_{\textnormal{eq}}(\boldsymbol \eta) \Lambda^1_{\mu,ijk}\, G^{-1}_{\mu,\nu}\, \Lambda^1_{\nu, mno} \sqrt g\, \textnormal d \boldsymbol\eta.\]

where \(n^0_{\textnormal{eq}}(\boldsymbol \eta)\) is an MHD equilibrium density (0-form).

property WMM#

create_weighted_mass(V_id: str, W_id: str, *, name: str | None = None, weights: list | str | None = None, assemble: bool = False, transposed: bool = False)[source]#

Weighted mass matrix \(V^\alpha_h \to V^\beta_h\) with given (matrix-valued) weight function \(W(\boldsymbol \eta)\):

\[\mathbb M_{(\mu, ijk), (\nu, mno)}(W) = \int \Lambda^\beta_{\mu, ijk}\, W_{\mu,\nu}(\boldsymbol \eta)\, \Lambda^\alpha_{\nu, mno} \, \textnormal d \boldsymbol\eta.\]

Here, \(\alpha \in \{0, 1, 2, 3, v\}\) indicates the domain and \(\beta \in \{0, 1, 2, 3, v\}\) indicates the co-domain of the operator.

Parameters:

V_id (str) – Specifier for the domain of the operator (‘H1’, ‘Hcurl’, ‘Hdiv’, ‘L2’ or ‘H1vec’).
W_id (str) – Specifier for the co-domain of the operator (‘H1’, ‘Hcurl’, ‘Hdiv’, ‘L2’ or ‘H1vec’).
name (str) – Name of the operator.
weights (None | str | list) –
Information about the weights/block structure of the operator. Four cases are possible:
1. str : for square block matrices (V=W), a symmetry can be set in order to accelerate the assembly process. Possible strings are symm (symmetric), asym (anti-symmetric) and diag (diagonal).
2. None : all blocks are allocated, disregarding zero-blocks or any symmetry.
3. 1D list : 1d list consisting of either a) strings or b) matrices (3x3 callables or 3x3 list) and can be mixed. Predefined names are G, Ginv, DFinv, sqrt_g. Access them using strings in the 1d list: weights=['<name>']. Possible choices for key-value pairs in weights are, at the moment: eq_mhd: MHDequilibrium. To access them, use for <name> the string eq_<method name>, where <method name> can be found in the just mentioned base classes for MHD equilibria. By default, all scalars are multiplied. For division of scalars use 1/<name>.
4. 2D list : 2d list with the same number of rows/columns as the number of components of the domain/codomain spaces. The entries can be either a) callables or b) xp.ndarrays representing the weights at the quadrature points. If an entry is zero or None, the corresponding block is set to None to accelerate the dot product.
assemble (bool) – Whether to assemble the weighted mass matrix, i.e. computes the integrals with assemble
transposed (bool) – Whether to assemble the transposed operator.

Returns:

out

Return type:

A WeightedMassOperator object.

class H1vecMassMatrix_density(derham, mass_ops, domain)[source]#

Bases: object

Wrapper around a Weighted mass operator from H1vec to H1vec whose weights are given by a 3 form

property massop#: The WeightedMassOperator

property inv#: The inverse WeightedMassOperator

update_weight(coeffs)[source]#: Update the weighted mass matrix operator

class struphy.feec.mass.WeightedMassOperatorsOldForTesting(derham: Derham, domain: Domain, matrix_free: bool = False, **weights)[source]#

Bases: object

Collection of pre-defined struphy.feec.mass.WeightedMassOperator.

Parameters:

derham (Derham) – Discrete de Rham sequence on the logical unit cube.
domain (Geometry) – Mapping from logical unit cube to physical domain and corresponding metric coefficients.
**weights (dict) – Objects to access callables that can serve as weight functions.
matrix_free (bool) – If set to true will not compute the matrix associated with the operator but directly compute the product when called

Notes

Possible choices for key-value pairs in weights are, at the moment:

eq_mhd: MHDequilibrium

property derham#: Discrete de Rham sequence on the logical unit cube.

property domain#: Mapping from the logical unit cube to the physical domain with corresponding metric coefficients.

property weights#: Dictionary of objects that provide access to callables that can serve as weight functions.

property selected_weight#: String identifying one key of “weigths”. This key is used when selecting weight functions.

G(e1, e2, e3)[source]#: Metric tensor callable.

Ginv(e1, e2, e3)[source]#: Inverse metric tensor callable.

sqrt_g(e1, e2, e3)[source]#: Jacobian determinant callable.

DFinv(e1, e2, e3)[source]#: Inverse Jacobian callable.

property M0#: Mass matrix

\[\mathbb M^0_{ijk, mno} = \int \Lambda^0_{ijk}\, \Lambda^0_{mno} \sqrt g\, \textnormal d \boldsymbol\eta.\]

property M1#: Mass matrix

\[\mathbb M^1_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^1_{\mu,ijk}\, G^{-1}\, \vec{\Lambda}^1_{\nu, mno} \sqrt g\, \textnormal d \boldsymbol\eta.\]

property M2#: Mass matrix

\[\mathbb M^2_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^2_{\mu,ijk}\, G\, \vec{\Lambda}^2_{\nu, mno} \frac{1}{\sqrt g}\, \textnormal d \boldsymbol\eta.\]

property M3#: Mass matrix

\[\mathbb M^3_{ijk, mno} = \int \Lambda^3_{ijk}\, \Lambda^3_{mno} \frac{1}{\sqrt g}\, \textnormal d \boldsymbol\eta.\]

property Mv#: Mass matrix

\[\mathbb M^v_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^v_{\mu,ijk}\, G\, \vec{\Lambda}^v_{\nu, mno} \sqrt g\, \textnormal d \boldsymbol\eta.\]

property M1n#

Mass matrix

\[\mathbb M^{1,n}_{(\mu,ijk), (\nu,mno)} = \int n^0_{\textnormal{eq}}(\boldsymbol \eta) \vec{\Lambda}^1_{\mu,ijk}\, G^{-1}\, \vec{\Lambda}^1_{\nu, mno} \sqrt g\, \textnormal d \boldsymbol\eta.\]

where \(n^0_{\textnormal{eq}}(\boldsymbol \eta)\) is an MHD equilibrium density (0-form).

property M2n#

Mass matrix

\[\mathbb M^{2,n}_{(\mu,ijk), (\nu,mno)} = \int n^0_{\textnormal{eq}}(\boldsymbol \eta) \vec{\Lambda}^2_{\mu,ijk}\, G\, \vec{\Lambda}^2_{\nu, mno} \frac{1}{\sqrt g}\, \textnormal d \boldsymbol\eta.\]

where \(n^0_{\textnormal{eq}}(\boldsymbol \eta)\) is an MHD equilibrium density (0-form).

property Mvn#

Mass matrix

\[\mathbb M^{v,n}_{(\mu,ijk), (\nu,mno)} = \int n^0_{\textnormal{eq}}(\boldsymbol \eta) \vec{\Lambda}^v_{\mu,ijk}\, G\, \vec{\Lambda}^v_{\nu, mno} \sqrt g\, \textnormal d \boldsymbol\eta.\]

where \(n^0_{\textnormal{eq}}(\boldsymbol \eta)\) is an MHD equilibrium density (0-form).

property M1ninv#

Mass matrix

\[\mathbb M^{1,\frac{1}{n}}_{(\mu,ijk), (\nu,mno)} = \int \frac{1}{n^0_{\textnormal{eq}}(\boldsymbol \eta)} \vec{\Lambda}^1_{\mu,ijk}\, G^{-1}\, \vec{\Lambda}^1_{\nu, mno} \sqrt g\, \textnormal d \boldsymbol\eta.\]

where \(n^0_{\textnormal{eq}}(\boldsymbol \eta)\) is an MHD equilibrium density (0-form).

property M1J#

Mass matrix

\[\mathbb M^{1,J}_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^1_{\mu,ijk}\, G^{-1}\, \mathcal R^J\, \vec{\Lambda}^2_{\nu, mno} \, \textnormal d \boldsymbol\eta.\]

with the rotation matrix

\[\mathcal R^J_{\alpha, \nu} := \epsilon_{\alpha \beta \nu}\, J^2_{\textnormal{eq}, \beta}\,,\qquad s.t. \qquad \mathcal R^J \vec v = \vec J^2_{\textnormal{eq}} \times \vec v\,,\]

where \(\epsilon_{\alpha \beta \nu}\) stands for the Levi-Civita tensor and \(J^2_{\textnormal{eq}, \beta}\) is the \(\beta\)-component of the MHD equilibrium current density (2-form).

property M2J#

Mass matrix

\[\mathbb M^{2,J}_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^2_{\mu,ijk}\, \mathcal R^J\, \vec{\Lambda}^2_{\nu, mno} \, \frac{1}{\sqrt g}\, \textnormal d \boldsymbol\eta.\]

with the rotation matrix

\[\mathcal R^J_{\alpha, \nu} := \epsilon_{\alpha \beta \nu}\, J^2_{\textnormal{eq}, \beta}\,,\qquad s.t. \qquad \mathcal R^J \vec v = \vec J^2_{\textnormal{eq}} \times \vec v\,,\]

where \(\epsilon_{\alpha \beta \nu}\) stands for the Levi-Civita tensor and \(J^2_{\textnormal{eq}, \beta}\) is the \(\beta\)-component of the MHD equilibrium current density (2-form).

property MvJ#

Mass matrix

\[\mathbb M^{v,J}_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^v_{\mu,ijk}\, \mathcal R^J\, \vec{\Lambda}^v_{\nu, mno} \, \textnormal d \boldsymbol\eta.\]

with the rotation matrix

\[\mathcal R^J_{\alpha, \nu} := \epsilon_{\alpha \beta \nu}\, J^2_{\textnormal{eq}, \beta}\,,\qquad s.t. \qquad \mathcal R^J \vec v = \vec J^2_{\textnormal{eq}} \times \vec v\,,\]

where \(\epsilon_{\alpha \beta \nu}\) stands for the Levi-Civita tensor and \(J^2_{\textnormal{eq}, \beta}\) is the \(\beta\)-component of the MHD equilibrium current density (2-form).

property M2B_div0#

Mass matrix

\[\mathbb M^{2,B}_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^2_{\mu,ijk}\, \mathcal R^J\, \vec{\Lambda}^2_{\nu, mno} \, \frac{1}{\sqrt g}\, \textnormal d \boldsymbol\eta.\]

with the rotation matrix

\[\mathcal R^J_{\alpha, \nu} := \epsilon_{\alpha \beta \nu}\, B^2_{\textnormal{eq}, \beta}\,,\qquad s.t. \qquad \mathcal R^J \vec v = \vec B^2_{\textnormal{eq}} \times \vec v\,,\]

where \(\epsilon_{\alpha \beta \nu}\) stands for the Levi-Civita tensor and \(B^2_{\textnormal{eq}, \beta}\) is the \(\beta\)-component of the MHD equilibrium magnetic field (2-form).

property M2B#

Mass matrix

\[\mathbb M^{2,B}_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^2_{\mu,ijk}\, \mathcal R^J\, \vec{\Lambda}^2_{\nu, mno} \, \frac{1}{\sqrt g}\, \textnormal d \boldsymbol\eta.\]

with the rotation matrix

\[\mathcal R^J_{\alpha, \nu} := \epsilon_{\alpha \beta \nu}\, B^2_{\textnormal{eq}, \beta}\,,\qquad s.t. \qquad \mathcal R^J \vec v = \vec B^2_{\textnormal{eq}} \times \vec v\,,\]

where \(\epsilon_{\alpha \beta \nu}\) stands for the Levi-Civita tensor and \(B^2_{\textnormal{eq}, \beta}\) is the \(\beta\)-component of the MHD equilibrium magnetic field (2-form).

property M2Bn#

Mass matrix

\[\mathbb M^{2,BN}_{(\mu,ijk), (\nu,mno)} = \int \vec{\Lambda}^2_{\mu,ijk}\, \mathcal R^J\, \vec{\Lambda}^2_{\nu, mno} \, \frac{1}{n^0_{\textnormal{eq}}(\boldsymbol \eta)}\, \frac{1}{\sqrt g}\, \textnormal d \boldsymbol\eta.\]

with the rotation matrix

\[\mathcal R^J_{\alpha, \nu} := \epsilon_{\alpha \beta \nu}\, B^2_{\textnormal{eq}, \beta}\,,\qquad s.t. \qquad \mathcal R^J \vec v = \vec B^2_{\textnormal{eq}} \times \vec v\,,\]

where \(\epsilon_{\alpha \beta \nu}\) stands for the Levi-Civita tensor and \(B^2_{\textnormal{eq}, \beta}\) is the \(\beta\)-component of the MHD equilibrium magnetic field (2-form).

property M1Bninv#

Mass matrix

\[\mathbb M^{1,B\frac{1}{n}}_{(\mu,ijk), (\nu,mno)} = \int \frac{1}{n^0_{\textnormal{eq}}(\boldsymbol \eta)}\, \vec{\Lambda}^1_{\mu,ijk}\, G^{-1}\, \mathcal R^J_{\alpha, \gamma}\, G^{-1}_{\gamma,\nu}\, \vec{\Lambda}^1_{\nu, mno} \, \sqrt g\, \textnormal d \boldsymbol\eta.\]

with the rotation matrix

\[\mathcal R^J_{\alpha, \nu} := \epsilon_{\alpha \beta \nu}\, B^2_{\textnormal{eq}, \beta}\,,\qquad s.t. \qquad \mathcal R^J \vec v = \vec B^2_{\textnormal{eq}} \times \vec v\,,\]

where \(\epsilon_{\alpha \beta \nu}\) stands for the Levi-Civita tensor and \(B^2_{\textnormal{eq}, \beta}\) is the \(\beta\)-component of the MHD equilibrium magnetic field (2-form).

property M1perp#: Mass matrix

\[\begin{split}\mathbb M^{1, \perp}_{(\mu, ijk), (\nu, mno)} = \int \vec{\Lambda}^1_{\mu, ijk}\, DF^{-1} \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix} DF^{-\top} \vec{\Lambda}^1_{\nu, mno} \sqrt g\, \textnormal d \boldsymbol\eta.\end{split}\]

property M0ad#: Mass matrix

\[\mathbb M^0_{ijk, mno} = \int \Lambda^0_{ijk}\, \Lambda^0_{mno} \sqrt g\, \textnormal d \boldsymbol\eta.\]

property M1gyro#

Mass matrix

\[\mathbb M^{1,n}_{(\mu,ijk), (\nu,mno)} = \int n^0_{\textnormal{eq}}(\boldsymbol \eta) \Lambda^1_{\mu,ijk}\, G^{-1}_{\mu,\nu}\, \Lambda^1_{\nu, mno} \sqrt g\, \textnormal d \boldsymbol\eta.\]

where \(n^0_{\textnormal{eq}}(\boldsymbol \eta)\) is an MHD equilibrium density (0-form).

Bases: LinOpWithTransp

Class for assembling weighted mass matrices in 3d.

Weighted mass matrices \(\mathbb M^{\beta\alpha}: \mathbb R^{N_\alpha} \to \mathbb R^{N_\beta}\) are of the general form

\[\mathbb M^{\beta \alpha}_{(\mu,ijk),(\nu,mno)} = \int_{[0, 1]^3} \Lambda^\beta_{\mu,ijk} \, A_{\mu,\nu} \, \Lambda^\alpha_{\nu,mno} \, \textnormal d^3 \boldsymbol\eta\,,\]

where the weight fuction \(A\) is a tensor of rank 0, 1 or 2, depending on domain and co-domain of the operator, and \(\Lambda^\alpha_{\nu, mno}\) is the B-spline basis function with tensor-product index \(mno\) of the \(\nu\)-th component in the space \(V^\alpha_h\). These matrices are sparse and stored in StencilMatrix format.

Finally, \(\mathbb M^{\beta\alpha}\) can be multiplied by PolarExtractionOperator and BoundaryOperator, \(\mathbb B\, \mathbb E\, \mathbb M^{\beta\alpha} \mathbb E^T \mathbb B^T\), to account for Polar splines and/or FEEC boundary conditions, respectively.

Parameters:

derham (Derham) – Struphy Derham object.
V (TensorFemSpace | VectorFemSpace) – Tensor product spline space from psydac.fem.tensor (domain, input space).
W (TensorFemSpace | VectorFemSpace) – Tensor product spline space from psydac.fem.tensor (codomain, output space).
name (str) – Name of the operator.
V_extraction_op (PolarExtractionOperator, optional) – Extraction operator to polar sub-space of V.
W_extraction_op (PolarExtractionOperator, optional) – Extraction operator to polar sub-space of W.
V_boundary_op (BoundaryOperator, optional) – Boundary operator that sets essential boundary conditions.
W_boundary_op (BoundaryOperator, optional) – Boundary operator that sets essential boundary conditions.
weights_info (NoneType | str | list) –
Information about the weights/block structure of the operator. Three cases are possible:
1. None : all blocks are allocated, disregarding zero-blocks or any symmetry.
2. str : for square block matrices (V=W), a symmetry can be set in order to accelerate the assembly process. Possible strings are symm (symmetric), asym (anti-symmetric) and diag (diagonal).
3. list : 2d list with the same number of rows/columns as the number of components of the domain/codomain spaces. The entries can be either a) callables or b) xp.ndarrays representing the weights at the quadrature points. If an entry is zero or None, the corresponding block is set to None to accelerate the dot product.
transposed (bool) – Whether to assemble the transposed operator.
matrix_free (bool) – If set to true will not compute the matrix associated with the operator but directly compute the product when called

property derham#

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 name#

property domain_femspace#

property codomain_femspace#

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 M#

property M0#

property matrix#

property nquads#

property symmetry#

property weights#

dot(v, out=None, apply_bc=True)[source]#

Dot product of the operator with a vector.

Parameters:

v (psydac.linalg.basic.Vector) – The input (domain) vector.
out (psydac.linalg.basic.Vector, optional) – If given, the output will be written in-place into this vector.
apply_bc (bool) – Whether to apply the boundary operators (True) or not (False).

Returns:

out – The output (codomain) vector.

Return type:

psydac.linalg.basic.Vector

transpose(conjugate=False)[source]#: Returns the transposed operator.

assemble(weights=None, clear=True, verbose=True)[source]#

Assembles the weighted mass matrix, i.e. computes the integrals

\[\mathbb M^{\beta \alpha}_{(\mu,ijk),(\nu,mno)} = \int_{[0, 1]^3} \Lambda^\beta_{\mu,ijk} \, A_{\mu,\nu} \, \Lambda^\alpha_{\nu,mno} \, \textnormal d^3 \boldsymbol\eta\,.\]

The integration is performed with Gauss-Legendre quadrature over the logical domain.

Parameters:

weights (list | NoneType) – Weight function(s) (callables or xp.ndarrays) in a 2d list of shape corresponding to number of components of domain/codomain. If weights=None, the weight is taken from the given weights in the instanziation of the object, else it will be overriden.
clear (bool) – Whether to first set all data to zero before assembly. If False, the new contributions are added to existing ones.
verbose (bool) – Whether to do some printing.

copy(out=None)[source]#

Create a copy of self, that can potentially be stored in a given WeightedMassOperator.

Parameters:: out (WeightedMassOperator(optional)) – The existing WeightedMassOperator in which we want to copy self.

static eval_quad(W, coeffs, out=None)[source]#

Evaluates a given FEM field defined by its coefficients at the L2 quadrature points.

Parameters:

W (TensorFemSpace | VectorFemSpace) – Tensor product spline space from psydac.fem.tensor.
coeffs (StencilVector | BlockVector) – The coefficient vector corresponding to the FEM field. Ghost regions must be up-to-date!
out (xp.ndarray | list/tuple of xp.ndarrays, optional) – If given, the result will be written into these arrays in-place. Number of outs must be compatible with number of components of FEM field.

Returns:

out – The values of the FEM field at the quadrature points.

Return type:

xp.ndarray | list/tuple of xp.ndarrays

class struphy.feec.mass.StencilMatrixFreeMassOperator(derham, V, W, weights=None, nquads=None)[source]#

Bases: LinOpWithTransp

Class implementing matrix-free weighted mass operators between StencilVectorSpaces.

The result of the dot product with a spline function \(S_h\) is computed as

\[w^\mu_{ijk} = \int \Lambda_{\mu,ijk}\, S_h\, w(\boldsymbol\eta)\,\textrm d \boldsymbol \eta \,,\]

where \(w(\boldsymbol\eta)\) is a weight function (including the geometric weights).

Should only be instanciated via WeightedMassOperator, where it’s used to replace StencilMatrix when one does not want to assemble the matrix for cost reasons

Parameters:

V (TensorFemSpace) – Domain of the mass operator
W (TensorFemSpace) – Codomain of the mass operator
weights (callable | numpy.ndarry | None) – The weights of the mass operator

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 nquads#

property derham#: Discrete de Rham sequence on the logical unit cube.

property tosparse#: Convert to a sparse matrix in any of the formats supported by scipy.sparse.

property toarray#: Convert to Numpy 2D array.

transpose(conjugate=False)[source]#

Transpose the LinearOperator .

If conjugate is True, return the Hermitian transpose.

property weights#

dot(v, out=None)[source]#

Dot product of the operator with a vector. Direct computation (not using a StencilMatrix).

Parameters:

v (psydac.linalg.basic.Vector) – The input (domain) vector.
out (psydac.linalg.basic.Vector, optional) – If given, the output will be written in-place into this vector.
apply_bc (bool) – Whether to apply the boundary operators (True) or not (False).

Returns:

out – The output (codomain) vector.

Return type:

psydac.linalg.basic.Vector

diagonal(inverse=False, sqrt=False, out=None)[source]#

Get the coefficients on the main diagonal as a StencilDiagonalMatrix object.

Parameters:

inverse (bool) – If True, get the inverse of the diagonal. (Default: False).
sqrt (bool) – If True, get the square root of the diagonal. (Default: False). Can be combined with inverse to get the inverse square root
out (StencilDiagonalMatrix) – If provided, write the diagonal entries into this matrix. (Default: None).

Returns:

The matrix which contains the main diagonal of self (or its inverse).

Return type:

StencilDiagonalMatrix

Table of Contents

Weighted mass operators#

This Page