Is there any built-in method in FEniCS legacy to calculate the inverse square root of FEniCS matrices?

Legacy DOLFIN
1

from fenics import *
from scipy.linalg import fractional_matrix_power

def compute_tau_supg(dt, A1, A2, K11, K12, K21, K22, mesh):
“”"
Compute τ_SUPG according to:
τ_SUPG = (4/Δt² I + GᵢⱼAᵢAⱼ + CₗGᵢⱼGₖₗKᵢₖKₗⱼ)^(-1/2)

Parameters:
    dt: time step size
    A1, A2: Advective flux Jacobian matrices
    K11, K12, K21, K22: Diffusion matrices
    mesh: FEniCS mesh
"""

# Get spatial coordinates and compute metric tensor Gᵢⱼ
x = SpatialCoordinate(mesh)
J = grad(x)  # Jacobian matrix
G = as_matrix([
    [dot(J[:, 0], J[:, 0]), dot(J[:, 0], J[:, 1])],
    [dot(J[:, 0], J[:, 1]), dot(J[:, 1], J[:, 1])]
])

# Create identity matrix of same size as advection matrices
I = Identity(A1.ufl_shape[0])

# Compute GᵢⱼAᵢAⱼ term
AA_term = (G[0, 0] * A1 * A1 +
           G[0, 1] * (A1 * A2 + A2 * A1) +
           G[1, 1] * A2 * A2)

# Compute CₗGᵢⱼGₖₗKᵢₖKₗⱼ term
C_I = Constant(1.0)  # Inverse estimate constant
KK_term = C_I * (
        G[0, 0] * G[0, 0] * (K11 * K11 + K12 * K12) +
        G[0, 1] * G[0, 1] * (K11 * K21 + K12 * K22) +
        G[1, 0] * G[1, 0] * (K21 * K11 + K22 * K12) +
        G[1, 1] * G[1, 1] * (K21 * K21 + K22 * K22)
)

# Compute matrix under square root: 4/Δt² I + GᵢⱼAᵢAⱼ + CₗGᵢⱼGₖₗKᵢₖKₗⱼ
M = 4.0 / (dt * dt) * I + AA_term + KK_term

# Compute τ_SUPG using matrix power
tau_SUPG = as_matrix(fractional_matrix_power(M, -0.5))

return tau_SUPG

this is code, I assemble the matrices used as parameters by using as_matrix, then I want to calculate the inverse square root of these matrices, is there any built-in function to do this?

2

Dear @Deep_Math,
This seems to be a continuation of How to compute the (-0.5) power of a matrix - #2 by Stein

where you haven’t answered the question provided by other users.

Secondly, you problem is not reproducible, as we don’t know the size of your matrix. Additionally, you haven’t provided the necessary definitions to execut you code.

Thirdly, I assume that there are ufl components within the matrix, such as spatial coordinates (which are symbolic entities, not numbers), and can thus not be inverted with scipy.

Note that for 2x2 and3x3 matrices there is a symbolical inverse in ufl, but you cannot take any power of a matrix in ufl (except 2).

What I think you could do in DOLFINx is to use the @a-latyshev external operator library (GitHub - a-latyshev/dolfinx-external-operator: Extension of DOLFINx implementing the concept of external operator ) to get a place to patch in the scipy fractional matrix implementation (which is based on: https://epubs.siam.org/doi/abs/10.1137/10081232X)