Efficient assembling of a matrix out of an existing matrix

Question

I am trying to implement a FEM-FCT algorithm for solving a transport equation for a special purpose. Source of the algorithm: http://www.mathematik.uni-dortmund.de/~kuzmin/linfct.pdf

My specific problem is, that I have a known Matrix $K_{ij}$ and I have to assemble a new matrix $D_{ij}$ using components of $K_{ij}$ the following way:

$$D_{ij}=max(-K_{ij},0,-K_{ji})$$

1. Is there a efficient/easy way to do this?

2. I have found a dirty way to achieve it, using the set command. But I guess it is slow and it also seems that there is a bug:
   I write blocks of components into the new matrix. It is working for a 3x3 block, while it produces an error for others, like 4x4. This seems to me like it is not working as intended? (See code below.)

from fenics import *
import numpy as np
import os

mesh=RectangleMesh(Point(0,0), Point(1, 1), 1, 1, "right/left")

dA = Measure("ds", domain=mesh)
dV = Measure("dx", domain=mesh)


ScalarSpace = FunctionSpace(mesh,'CG',1)
VectorSpace = VectorFunctionSpace(mesh,'CG',1)

test_rho = TestFunction(ScalarSpace)
delta_rho = TrialFunction(ScalarSpace)

v_g_exp = Expression(('1.0','0.0'),degree=2, domain=mesh)
v_g = interpolate(v_g_exp, VectorSpace)


####################  K and D matrix ###########################

k_form = div(delta_rho*v_g)*test_rho*dV
K = assemble(k_form)
D = assemble(k_form)
D.zero()
print K.array()

####################  working  #################################
print D
block = np.ones(9,dtype=np.float_)
print block
rows = np.array([0,1,2],dtype=np.intc)
print rows
cols = np.array([0,1,2],dtype=np.intc)
print cols
D.set(block,rows,cols)
D.apply('add')
print D.array()

####################  not working  ############################
print D
block = np.ones(16,dtype=np.float_)
print block
rows = np.array([0,1,2,3],dtype=np.intc)
print rows
cols = np.array([0,1,2,3],dtype=np.intc)
print cols
D.set(block,rows,cols)
D.apply('add')
print D.array()

###############################################################

(3. How will this react regarding a potential future parallelization?)

Similar question was asked here: https://fenicsproject.org/qa/9017/fem-fct-artificial-diffusion-matrix-construction

Answer 1

The easiest way to do this is using backend functions for accessing matrix data. In particular, you want to avoid looping over matrix entries in the Python code.

The following example shows how you can do this if you have built dolfin with petsc4py.

K = assemble(k_form)

kmat = as_backend_type(K).mat()
dmat = kmat.duplicate()

# copy transpose of kmat into dmat
kmat.transpose(dmat)

# get values from kmat and dmat
i, c, kvals = kmat.getValuesCSR()
_, _, dvals = dmat.getValuesCSR()

# compute values for matrix D
dvals = - reduce(numpy.minimum, (kvals, dvals, 0))

# copy values into matrix D
dmat.setValuesCSR(i, c, dvals)
dmat.assemble()
D = PETScMatrix(dmat)

Comments

Helped a lot, thank you.