Preconditioners

class struphy.feec.preconditioner.MassMatrixPreconditioner(mass_operator, apply_bc=True, dim_reduce=0)

Bases: LinearOperator

Preconditioner for inverting 3d weighted mass matrices.

The mass matrix is approximated by a Kronecker product of 1d mass matrices in each direction with correct boundary conditions (block diagonal in case of vector-valued spaces). In this process, the 3d weight function is appoximated by a 1d counterpart in the dim_reduce direction (default 1st direction) at the fixed point (0.5) in the other directions. The inversion is then performed with a Kronecker solver.

Parameters:

• mass_operator (struphy.feec.mass.WeightedMassOperator) – The weighted mass operator for which the approximate inverse is needed.

• apply_bc (bool) – Whether to include boundary operators.

• dim_reduce (int) – Along which axis to take the approximate value of the weight

property space

Stencil-/BlockVectorSpace or PolarDerhamSpace.

property matrix

Approximation of the input mass matrix as KroneckerStencilMatrix.

property solver

KroneckerLinearSolver or BlockDiagonalSolver for exactly inverting the approximate mass matrix self.matrix.

property codomain

The codomain of the linear operator - an element of Vectorspace

property domain

The domain of the linear operator - an element of Vectorspace

property dtype

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

tosparse()

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

toarray()

Convert to Numpy 2D array.

transpose(conjugate=False)

Returns the transposed operator.

solve(rhs, out=None)

Computes (B * E * M^(-1) * E^T * B^T) * rhs as an approximation for an inverse mass matrix.

Parameters:

• rhs (psydac.linalg.basic.Vector) – The right-hand side vector.

• out (psydac.linalg.basic.Vector, optional) – If given, the output vector will be written into this vector in-place.

Returns:

out – The result of (B * E * M^(-1) * E^T * B^T) * rhs.

Return type:

psydac.linalg.basic.Vector

dot(v, out=None)

Apply linear operator to Vector v. Result is written to Vector out, if provided.

class struphy.feec.preconditioner.MassMatrixDiagonalPreconditioner(mass_operator, apply_bc=True)

Bases: LinearOperator

Preconditioner for inverting 3d weighted mass matrices. The mass matrix is approximated by

\[D^{1/2} * \hat D^{-1/2} * \hat M * \hat D^{-1/2} * D^{1/2}\]

Where $D$ is the diagonal of the matrix to invert, \(\hat M\) is the mass matrix on the logical domain that is a Kronecker product (fastly inverted) and \(\hat D^{-1/2}\) is the diagonal of \(\hat M\).

Notes

Reference: G. Loli, G. Sangalli, M. Tani, “Easy and efficient preconditioning of the isogeometric mass matrix”, Comp. Math. Appl., Vol. 116, 2022

Parameters:

• mass_operator (WeightedMassOperator) – The weighted mass operator for which the approximate inverse is needed.

• apply_bc (bool) – Whether to include boundary operators.

property space

Stencil-/BlockVectorSpace or PolarDerhamSpace.

property matrix

Mass matrix on the logical domain as KroneckerStencilMatrix.

property solver

KroneckerLinearSolver or BlockDiagonalSolver for exactly inverting the approximate mass matrix self.matrix.

property codomain

The codomain of the linear operator - an element of Vectorspace

property domain

The domain of the linear operator - an element of Vectorspace

property dtype

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

update_mass_operator(mass_operator)

Update the mass operator to enable recycling the preconditioner

tosparse()

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

toarray()

Convert to Numpy 2D array.

transpose(conjugate=False)

Returns the transposed operator.

solve(rhs, out=None)

Computes \((B * E * M^{-1} * E^T * B^T) * rhs\) as an approximation for an inverse mass matrix, with \(M = D^{1/2} * \hat D^{-1/2} * \hat M * \hat D^{-1/2} * D^{1/2}\).

Parameters:

• rhs (Vector) – The right-hand side vector.

• out (Vector, optional) – If given, the output vector will be written into this vector in-place.

Returns:

out – The result of \((B * E * M^{-1} * E^T * B^T) * rhs\).

Return type:

Vector

dot(v, out=None)

Apply linear operator to Vector v. Result is written to Vector out, if provided.

class struphy.feec.preconditioner.FFTSolver(circmat)

Bases: BandedSolver

Solve the equation Ax = b for x, assuming A is a circulant matrix. b can contain multiple right-hand sides (RHS) and is of shape (#RHS, N).

Parameters:

circmat (xp.ndarray) – Generic circulant matrix.

property space

solve(rhs, out=None, transposed=False)

Solves for the given right-hand side.

Parameters:

• rhs (xp.ndarray) – The right-hand sides to solve for. The vectors are assumed to be given in C-contiguous order, i.e. if multiple right-hand sides are given, then rhs is a two-dimensional array with the 0-th index denoting the number of the right-hand side, and the 1-st index denoting the element inside a right-hand side.

• out (xp.ndarray, optional) – Output vector. If given, it has to have the same shape and datatype as rhs.

• transposed (bool) – If and only if set to true, we solve against the transposed matrix. (supported by the underlying solver)

struphy.feec.preconditioner.is_circulant(mat)

Returns true if a matrix is circulant.

Parameters:

mat (array[float]) – The matrix that is checked to be circulant.

Returns:

circulant – Whether the matrix is circulant (=True) or not (=False).

Return type:

bool

Linear operators

class psydac.linalg.basic.LinearOperator

Bases: ABC

Abstract base class for all linear operators acting between two vector spaces V (domain) and W (codomain).

property shape

A tuple containing the dimension of the codomain and domain.

abstract property domain

The domain of the linear operator - an element of Vectorspace

abstract property codomain

The codomain of the linear operator - an element of Vectorspace

abstract property dtype

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

abstract tosparse()

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

abstract toarray()

Convert to Numpy 2D array.

abstract dot(v, out=None)

Apply the LinearOperator self to the Vector v.

The result is written to the Vector out, if provided.

Parameters:

• v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.

• out (Vector) – The vector in which the result of the operation is stored. It must belong to the codomain of self. If out is None, a new vector is created and returned.

Returns:

The result of the operation. If out is None, a new vector is returned. Otherwise, the result is stored in out and out is returned.

Return type:

Vector

abstract transpose(conjugate=False)

Transpose the LinearOperator .

If conjugate is True, return the Hermitian transpose.

property T

Calls transpose method to return the transpose of self.

property H

Calls transpose method with conjugate=True flag to return the Hermitian transpose of self.

idot(v, out)

Implements out += self @ v with a temporary. Subclasses should provide an implementation without a temporary.

dot_inner(v, w)

Compute the inner product of (self @ v) with w, without a temporary.

This is equivalent to self.dot(v).inner(w), but avoids the creation of a temporary vector because the result of self.dot(v) is stored in a local work array. If self is a positive-definite operator, this operation is a (weighted) inner product.

Parameters:

• v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.

• w (Vector) – The second vector in the inner product. It must belong to the codomain of self.

Returns:

The result of the inner product between (self @ v) and w. If the field of self is real, this is a real number. If the field of self is complex, this is a complex number.

Return type:

float | complex

class psydac.linalg.block.BlockLinearOperator(V1, V2, blocks=None)

Bases: LinearOperator

Linear operator that can be written as blocks of other Linear Operators. Either the domain or the codomain of this operator, or both, should be of class BlockVectorSpace.

Parameters:

• V1 (psydac.linalg.block.VectorSpace) – Domain of the new linear operator.

• V2 (psydac.linalg.block.VectorSpace) – Codomain of the new linear operator.

• blocks (dict | (list of lists) | (tuple of tuples)) –

  LinearOperator objects (optional).

  1. ’blocks’ can be dictionary with

     . key = tuple (i, j), where i and j are two integers >= 0 . value = corresponding LinearOperator Lij

  2. ’blocks’ can be list of lists (or tuple of tuples) where blocks[i][j]

     is the LinearOperator Lij (if None, we assume null 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.

conjugate(out=None)

conj(out=None)

tosparse(**kwargs)

Convert to any Scipy sparse matrix format.

toarray(**kwargs)

Convert to Numpy 2D array.

dot(v, out=None)

Apply the LinearOperator self to the Vector v.

The result is written to the Vector out, if provided.

Parameters:

• v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.

• out (Vector) – The vector in which the result of the operation is stored. It must belong to the codomain of self. If out is None, a new vector is created and returned.

Returns:

The result of the operation. If out is None, a new vector is returned. Otherwise, the result is stored in out and out is returned.

Return type:

Vector

transpose(conjugate=False, out=None)

” Return the transposed BlockLinearOperator, or the Hermitian Transpose if conjugate==True

Parameters:

• conjugate (Bool(optional)) – True to get the Hermitian adjoint.

• out (BlockLinearOperator(optional)) – Optional out for the transpose to avoid temporaries

diagonal(*, inverse=False, sqrt=False, out=None)

Get the coefficients on the main diagonal as another BlockLinearOperator object.

Parameters:

• inverse (bool) – If True, get the inverse of the diagonal. (Default: False). Can be combined with sqrt to get the inverse square root.

• 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 (BlockLinearOperator) – 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:

BlockLinearOperator

property blocks

Immutable 2D view (tuple of tuples) of the linear operator, including the empty blocks as ‘None’ objects.

property n_block_rows

property n_block_cols

property nonzero_block_indices

Tuple of (i, j) pairs which identify the non-zero blocks: i is the row index, j is the column index.

update_ghost_regions()

exchange_assembly_data()

remove_spurious_entries()

property ghost_regions_in_sync

transform(operation)

Applies an operation on each block in this BlockLinearOperator.

Parameters:

operation (LinearOperator -> LinearOperator) – The operation which transforms each block.

backend()

copy(out=None)

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

Parameters:

out (BlockLinearOperator(optional)) – The existing BlockLinearOperator in which we want to copy self.

Returns:

The copy of self, either stored in the given BlockLinearOperator out (if provided) or in a new one. In the corner case where out=self the self object is immediately returned.

Return type:

BlockLinearOperator

topetsc()

Convert to petsc data structure.

compute_interface_matrices_transpose()

set_backend(backend, precompiled=False)

class struphy.feec.linear_operators.LinOpWithTransp

Bases: LinearOperator

Base class for linear operators that MUST implement a transpose method which returns a new transposed operator.

abstract transpose()

Transpose the LinearOperator .

If conjugate is True, return the Hermitian transpose.

toarray_struphy(out=None, is_sparse=False, format='csr')

Transforms the linear operator into a matrix, which is either stored in dense or sparse format.

Parameters:

• out (Numpy.ndarray, optional) – If given, the output will be written in-place into this array.

• is_sparse (bool, optional) – If set to True the method returns the matrix as a Scipy sparse matrix, if set to false it returns the full matrix as a Numpy.ndarray

• format (string, optional) – Only relevant if is_sparse is True. Specifies the format in which the sparse matrix is to be stored. Choose from “csr” (Compressed Sparse Row, default),”csc” (Compressed Sparse Column), “bsr” (Block Sparse Row ), “lil” (List of Lists), “dok” (Dictionary of Keys), “coo” (COOrdinate format) and “dia” (DIAgonal).

Returns:

out – The matrix form of the linear operator. If ran in parallel each rank gets the full matrix representation of the linear operator.

Return type:

Numpy.ndarray or scipy.sparse.csr.csr_matrix

class struphy.feec.linear_operators.BoundaryOperator(vector_space, space_id, dirichlet_bc)

Bases: LinOpWithTransp

Applies homogeneous Dirichlet boundary conditions to a vector.

Parameters:

• vector_space (psydac.linalg.basic.VectorSpace) – The vector space associated to the operator.

• space_id (str) – Symbolic space ID of vector_space (H1, Hcurl, Hdiv, L2 or H1vec).

• dirichlet_bc (tuple[tuple[bool]]) – Whether to apply homogeneous Dirichlet boundary conditions (at left or right boundary in each direction).

property domain

The domain of the linear operator - an element of Vectorspace

property codomain

The codomain of the linear operator - an element of Vectorspace

property dtype

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

property tosparse

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

property toarray

Convert to Numpy 2D array.

property bc

property dim_nz_pol

property dim_nz_tor

property dim_nz

dot(v, out=None)

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.

Returns:

out – The output (codomain) vector.

Return type:

psydac.linalg.basic.Vector

transpose(conjugate=False)

Returns the transposed operator.

Utilities

class struphy.feec.utilities.RotationMatrix(*vec_fun)

Bases: object

For a given vector-valued function a(e1, e2, e3), creates the callable matrix R(e1, e2, e3) that represents the local rotation Rv = a x v at (e1, e2, e3) for any vector v in R^3.

When called, R(e1, e2, e3) is a five-dimensional array, with the 3x3 matrix in the last two indices.

Parameters:

*vec_fun (list) – Three callables that represent the vector-valued function a.

struphy.feec.utilities.create_equal_random_arrays(V, seed=123, flattened=False)

Creates two equal random arrays, where one array is a numpy array and the other one a distributed psydac array.

Parameters:

• V (TensorFemSpace or VectorFemSpace) – The FEM space to which the arrays belong to.

• seed (int) – The seed used in the radom number generator.

• flattened (bool) – Whether to return flattened arrays or 3d arrays.

Returns:

• arr (list or array (if flattened)) – The 3d numpy arrays for each component of the space.

• arr_psy (StencilVector of BlockVector) – The distributed psydac array.

struphy.feec.utilities.compare_arrays(arr_psy, arr, rank, atol=1e-14, verbose=False)

Assert equality of distributed psydac array and corresponding fraction of cloned numpy array. Arrays can be block-structured as nested lists/tuples.

Parameters:

• arr_psy (psydac object) – Stencil/Block Vector/Matrix.

• arr (array) – Numpy array with same tuple/list structure as arr_psy. If arr is a matrix it can be stencil or band format.

• rank (int) – Rank of mpi process

• atol (float) – Absolute tolerance used in numpy.allclose()

struphy.feec.utilities.apply_essential_bc_to_array(space_id: str, vector: Vector, bc: tuple)

Sets entries corresponding to boundary B-splines to zero.

Parameters:

• space_id (str) – The name of the continuous functions space the given vector belongs to (H1, Hcurl, Hdiv, L2 or H1vec).

• vector (Vector) – The vector whose boundary values shall be set to zero.

• bc (tuple[tuple[bool]]) – Whether to apply homogeneous Dirichlet boundary conditions (at left or right boundary in each direction).

struphy.feec.utilities.create_weight_weightedmatrix_hybrid(b, weight_pre, derham, accum_density, domain)

Creates weights needed for asembling matrix of hybrid model with kinetic ions and massless electrons

Parameters:

• b (Stencilvector) – finite element coefficients of magnetic fields.

• self._weight_pre (list of 3D arrays storing weights)

• derham (derham class)

• accum_density (StencilMatrix) – the density obtained by deposition of particles

• domain (domain class)

Returns:

self._weight_pre

Return type:

list of 3D arrays storing weights

Variational utilities

class struphy.feec.variational_utilities.BracketOperator(derham: Derham, u: BlockVector)

Bases: LinOpWithTransp

The linear map \(\mathbb R^{3N_0} \to \mathbb R^{3N_0}\),

\[\mathbf v \in \mathbb R^{3N_0} \mapsto \mathbf w = (w_{\mu,ijk})_{\mu,ijk} \in \mathbb R^{3N_0}\,,\]

defined by

\[w_{\mu,ijk} = \int \hat{\mathbf m}(\boldsymbol \eta)\, G\, [\mathbf v^\top \vec{\boldsymbol \Lambda}^v, \vec{\Lambda}^v_{\mu,ijk}] \,\sqrt g\, \textnormal d\boldsymbol \eta\,,\]

where \(\hat{\mathbf m}(\boldsymbol \eta)\) is a given vector-field, and with the usual vector-field bracket

\[[\mathbf v^\top \vec{\boldsymbol \Lambda}^v, \vec{\Lambda}^v_{\mu,ijk}] = \mathbf v^\top \vec{\boldsymbol \Lambda}^v \cdot \nabla \vec{\Lambda}^v_{\mu,ijk} - \vec{\Lambda}^v_{\mu,ijk} \cdot \nabla (\mathbf v^\top \vec{\boldsymbol \Lambda}^v)\,.\]

This is discretized as

\[\mathbf w = \sum_{\mu = 1}^3 I_\mu \Big(\hat{\Pi}^{0}[\hat{\mathbf v}_h \cdot \vec{\boldsymbol \Lambda}^1 ] \mathbb G P_\mu - \hat{\Pi}^0[\hat{\mathbf A}^1_{\mu,h} \cdot \vec{\boldsymbol \Lambda}^v] \Big)^\top \mathbf u \,,\]

where \(I_\mu\) and \(P_\mu\) stand for the CoordinateInclusion and CoordinateProjector, respectively, and the vector \(\mathbf u = (\hat{\mathbf m}, \vec{\boldsymbol \Lambda}^v)_{L^2} = \mathbb M^v \mathbf m\) is provided as input. The weights in the the two BasisProjectionOperator are given by

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

Initialized and used in VariationalMomentumAdvection propagator.

Parameters:

• derham (Derham) – Discrete de Rham sequence.

• u (BlockVector) – Coefficient of a field belonging to the H1vec space of the de Rahm sequence, representing the mass matrix applie to the m factor in the above integral.

property domain

The domain of the linear operator - an element of Vectorspace

property codomain

The codomain of the linear operator - an element of Vectorspace

property dtype

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

property tosparse

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

property toarray

Convert to Numpy 2D array.

update_u(newu)

transpose(conjugate=False)

Transpose the LinearOperator .

If conjugate is True, return the Hermitian transpose.

dot(v, out=None)

Apply the LinearOperator self to the Vector v.

The result is written to the Vector out, if provided.

Parameters:

• v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.

• out (Vector) – The vector in which the result of the operation is stored. It must belong to the codomain of self. If out is None, a new vector is created and returned.

Returns:

The result of the operation. If out is None, a new vector is returned. Otherwise, the result is stored in out and out is returned.

Return type:

Vector

class struphy.feec.variational_utilities.L2_transport_operator(derham, transposed=False, weights=None)

Bases: LinOpWithTransp

Operator

\[\mathbf u \mapsto \nabla \cdot(\Pi^2(\rho \mathbf u)) \,\]

from H1vec to L2, where \(\rho\) is a discrete 3-form which can be updated.

Parameters:

• derham (Derham) – Discrete de Rham sequence.

• transposed (Bool) – Assemble the transposed 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 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)

Transpose the LinearOperator .

If conjugate is True, return the Hermitian transpose.

dot(v, out=None)

Apply the LinearOperator self to the Vector v.

The result is written to the Vector out, if provided.

Parameters:

• v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.

• out (Vector) – The vector in which the result of the operation is stored. It must belong to the codomain of self. If out is None, a new vector is created and returned.

Returns:

The result of the operation. If out is None, a new vector is returned. Otherwise, the result is stored in out and out is returned.

Return type:

Vector

update_coeffs(coeff)

Update the coefficient of the projection operator.

Parameters:

coeffs (StencilVector) – coefficient of the discrete 3 form to update the projection operator.

class struphy.feec.variational_utilities.Hdiv0_transport_operator(derham, transposed=False, weights=None)

Bases: LinOpWithTransp

Operator

\[u \mapsto \nabla \times (\Pi^1(\mathbf B \times \mathbf u)) \,\]

from H1vec to H(div), where \(\mathbf B\) is a discrete 2-form which can be updated.

Parameters:

• derham (Derham) – Discrete de Rham sequence.

• transposed (Bool) – Assemble the transposed 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 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)

Transpose the LinearOperator .

If conjugate is True, return the Hermitian transpose.

dot(v, out=None)

Apply the LinearOperator self to the Vector v.

The result is written to the Vector out, if provided.

Parameters:

• v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.

• out (Vector) – The vector in which the result of the operation is stored. It must belong to the codomain of self. If out is None, a new vector is created and returned.

Returns:

The result of the operation. If out is None, a new vector is returned. Otherwise, the result is stored in out and out is returned.

Return type:

Vector

update_coeffs(coeff)

Update the coefficient of the projection operator.

Parameters:

coeffs (BlockVector) – coefficient of the discrete 2 form to update the projection operator.

class struphy.feec.variational_utilities.Pressure_transport_operator(derham, phys_domain, Uv, gamma, transposed=False, weights1=None, weights2=None)

Bases: LinOpWithTransp

Operator

\[\mathbf u \mapsto \nabla \cdot (\Pi^2(p \mathbf u)) + (\gamma -1) \Pi^3(p \nabla \cdot \Pi^2(\mathbf u)) \,\]

from H1vec to L2, where \(p\) is a discrete 3-form which can be updated.

Parameters:

• derham (Derham) – Discrete de Rham sequence.

• phys_domain (Domain) – The domain in which the problem is discretized (needed for metric terms)

• Uv (BasisProjectionOperator) – The projection from H1vec to H(div)

• gamma (Float) – Thermodynamical constant

• transposed (Bool) – Assemble the transposed 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 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)

Transpose the LinearOperator .

If conjugate is True, return the Hermitian transpose.

dot(v, out=None)

Apply the LinearOperator self to the Vector v.

The result is written to the Vector out, if provided.

Parameters:

• v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.

• out (Vector) – The vector in which the result of the operation is stored. It must belong to the codomain of self. If out is None, a new vector is created and returned.

Returns:

The result of the operation. If out is None, a new vector is returned. Otherwise, the result is stored in out and out is returned.

Return type:

Vector

update_coeffs(coeff)

Update the coefficient of the projection operator.

Parameters:

coeffs (StencilVector) – coefficient of the discrete 3 form to update the projection operator.

class struphy.feec.variational_utilities.InternalEnergyEvaluator(derham, gamma)

Bases: object

Helper class for the evaluation of the internal energy or its partial derivative/discrete partial derivatives on an integration grid

This class only contains a lot of array corresponding to the integration grid to avoid the allocation of temporaries, and method that can be called to evaluate the energy and derivatives on the grid.

Parameters:

• derham (Derham) – Discrete de Rham sequence.

• gamma (Float) – Thermodynamical constant

ener(rho, s, out=None)

Themodynamical energy as a function of rho and s, usign the perfect gaz hypothesis.

\[E(\rho, s) = \rho^\gamma \text{exp}(s/\rho) \,.\]

dener_drho(rho, s, out=None)

Derivative with respect to rho of the thermodynamical energy as a function of rho and s, usign the perfect gaz hypothesis.

\[\frac{\partial E}{\partial \rho}(\rho, s) = (\gamma \rho^{\gamma-1} - s \rho^{\gamma-2})*\text{exp}(s/\rho) \,.\]

dener_ds(rho, s, out=None)

Derivative with respect to s of the thermodynamical energy as a function of rho and s, usign the perfect gaz hypothesis.

\[\frac{\partial E}{\partial s}(\rho, s) = \rho^{\gamma-1} \text{exp}(s/\rho) \,.\]

d2ener_drho2(rho, s, out=None)

Second derivative with respect to (rho, rho) of the thermodynamical energy as a function of rho and s, usign the perfect gaz hypothesis.

\[\frac{\partial^2 E}{\partial \rho^2}(\rho, s) = (\gamma*(\gamma-1) \rho^{\gamma-2}- 2 s (\gamma-1) rho^{\gamma-3}+ s^2 \rho^{\gamma-4}) \text{exp}(s/\rho) \,.\]

d2ener_ds2(rho, s, out=None)

Second derivative with respect to (s, s) of the thermodynamical energy as a function of rho and s, usign the perfect gaz hypothesis.

\[\frac{\partial^2 E}{\partial s^2}(\rho, s) = \rho^{\gamma-2} \text{exp}(s/ \rho) \,.\]

eta(delta_x, out=None)

Switch function \(\eta(\delta) = 1- \text{exp}((-\delta/10^{-5})^2)\).

evaluate_discrete_de_drho_grid(rhon, rhon1, sn, out=None)

Evaluate the discrete gradient of the internal energy with respect to the \(\rho\) variable

\[\eta(\delta \rho)\frac{e(\rho^{n+1},s^n)-e(\rho^{n},s^n)}{\rho^{n+1}-\rho^n}+(1-\eta(\delta \rho))\frac{\partial e}{\partial \rho}(\rho^{n+\frac{1}{2}}, s^n) \,,\]

evaluate_exact_de_drho_grid(rhon, sn, out=None)

Evaluation of the derivative of \(E\) with respect to \(\rho\) on the grid.

evaluate_discrete_de_ds_grid(rhon, sn, sn1, out=None)

Evaluate the discrete gradient of the internal energy with respect to the \(s\) variable

\[\eta(\delta \rho)\frac{e(\rho^{n},s^{n+1})-e(\rho^{n},s^n)}{s^{n+1}-s^n}+(1-\eta(\delta s))\frac{\partial e}{\partial s}(\rho^n, s^{n+\frac{1}{2}}) \,,\]

evaluate_exact_de_ds_grid(rhon, sn, out=None)

Evaluation of the derivative of \(E\) with respect to \(s\) on the grid.

evaluate_discrete_d2e_drho2_grid(rhon, rhon1, sn, out=None)

Evaluate the derivative of the discrete derivative with respect to rhon1

evaluate_discrete_d2e_ds2_grid(rhon, sn, sn1, out=None)

Evaluate the derivative of the discrete derivative with respect to sn1

eval_3form(coeffs, out=None)

Evaluate the 3 form with FE coefficient coeffs on the grid

class struphy.feec.variational_utilities.H1vecMassMatrix_density(derham, mass_ops, domain)

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)

Update the weighted mass matrix operator

class struphy.feec.variational_utilities.KineticEnergyEvaluator(derham, domain, mass_ops)

Bases: object

Helper class to evaluate the different Kinetic energy terms appearing in VariationalDensityEvolve.

This class only contains arrays corresponding to the integration grid to avoid the allocation of temporaries, methods that can be called to evaluate the energy and derivatives on the grid and weighted mass operators corresponding to integration against a vector field.

Parameters:

• derham (Derham) – Discrete de Rham sequence.

• domain (Domain) – The domain in which the problem is discretized (needed for metric terms)

• mass_ops (WeightedMassOperators) – The weighted mass operators needed to create new mass matrices

property M_un

Weighted mass matrix with domain H1vec et codomain L2 represented the integration against a vector field in H1vec

property M_un1

Weighted mass matrix with domain L2 et codomain H1vec represented the integration against a vector field in H1vec

get_u2_grid(un, un1, out)

Values of \(u_n \cdot u_{n+1}\) represented by the coefficient un and un1, on the integration grid

assemble_M_un(un)

Update the weights of the matrix M_un with the vector fields given by the coeficient un

assemble_M_un1(un1)

Update the weights of the matrix M_un1 with the vector fields given by the coeficient un1