Using ExternalOperator to Use Scipy Tools On Firedrake Matrix

[maggiemcc02]

Hello!

I define a matrix of functions (DG0) in my code but I then want to compute the square root of that matrix using scipy.linalg.sqrtm using an ExternalOperator. I am having trouble constructing an array I can give to scipy. Below is a general outline for the code. Any advice is appreciated!

class AbsoluteValueH(AbstractExternalOperator):

  def __init__(self, *operands, function_space, **kwargs):
      AbstractExternalOperator.__init__(self, *operands, function_space=function_space, **kwargs) 

  @assemble_method(0, (0,))
  def assemble_N(self, *args, **kwargs):
      """Evaluate N"""

      # my matrix : as_matrix([[func1, func2], [func3, func4]]) where funci (i=1,2,3,4) is in DG0 
      # of type list tensor 
      
      H = self.ufl_operands[0] 


      # Construct Hm that can be passed to scipy

      ... 

      # take the square root

      sqrtH = scipy.linalg.sqrtm(Hm)


      return sqrtH 

I apologize if this is a stupid question - I am still new-ish to Firedrake!

[connorjward]

If you have an assembled Firedrake Matrix then you can get an array representation of it if you run matrix.M.values. Note that that will not work in parallel and it returns a dense array so is not recommended for large problems.