I wish to build a sparse PETScMatrix (C, below) whose values are the point-wise values of a given function (rho). However I am unable to assign entries either element wise or through diagonal entries.
I cannot build a standard array of the values because at deployment mesh resolution a dense matrix contains too many values and causes memory problems. What would be the best way to accomplish this?
Is there a way to get directly a matrix representation C_ij = rho(x[i], x[j])
similar to the assemble command on a form?
from fenics import *
from mshr import *
import numpy as np
res = 20
r1 = 1
r0 = 0.8
origin = Point(0, 0)
domain = Circle(origin, r1) - Circle(origin, r0)
mesh = generate_mesh(domain, res)
V = FunctionSpace(mesh, "DG", 0)
tol = 1e-3
def rho(x):
val = np.exp(-5 * (np.power(x[0] - y[0], 2) + np.power(x[1] - y[1], 2)) )
if val < tol:
val = 0
return val
xs = V.tabulate_dof_coordinates().reshape(V.dim(), mesh.geometric_dimension())
C = PETScMatrix()
for i in range(len(xs)):
for j in range(len(xs)):
x = xs[i, :]
y = xs[j, :]
val = rho(x, y)
C.set(np.ascontiguousarray([val]), [i], [j])
Is there scope for generating C within a matrix-free algorithm? I may not understand your problem correctly but I can’t see how C can be anything but a dense matrix.
I am not sure if there is. Ultimately we want to be able to store the eigen-decomposition of C, making the random field it represents cheap to sample from. Would a matrix free representation allow us to build the eigen-expansion and save this to a file for access later?
That’s exactly what we’re trying to do! We want the resulting expansion to be quick to save / load since we use the random samples as a parameter in several simulations which we run in parallel.
It is possible to set up a scipy matrix and use it to create the PETScMatrix object for FEniCS. I modified your MWE, but the way you set that up definitely results in a dense matrix.
from dolfin import *
from mshr import *
import scipy.sparse as sp
import numpy as np
from petsc4py import PETSc
mesh = generate_mesh(Circle(Point(0, 0), 1) - Circle(Point(0, 0), 0.8), 20)
V = FunctionSpace(mesh, "DG", 0)
def rho(x, y):
val = np.exp(-5 * (np.power(x[0] - y[0], 2) + np.power(x[1] - y[1], 2)) )
if val < 1e-3:
val = 0
return val
xs = V.tabulate_dof_coordinates()
data = []
row_indices = []
col_indices = []
for i in range(len(xs)):
for j in range(len(xs)):
row_indices.append(i)
col_indices.append(j)
data.append(rho(xs[i, :], xs[j, :]))
C_intermediate = sp.csr_matrix( (data, (row_indices, col_indices)), shape=(V.dim(), V.dim()) )
C = PETScMatrix( PETSc.Mat().createAIJ(size = C_intermediate.shape, csr = (C_intermediate.indptr, C_intermediate.indices, C_intermediate.data)) )
Isn’t it a bad idea to use a sparse matrix data structure for a dense matrix though? What @jrtroy wants to compute isn’t a simple problem. The system is always dense. A matrix of around 3000 rows will probably hit the memory limit with a modern desktop. This means only very coarse approximations are possible using standard methods. There are many attempts to alleviate the computational costs, e.g., circulent embedding, matrix compression, H-matrix formulation or matrix-free assembly.
Just to add: If your domain is always going to be that annulus shape (two concentric circles). You may be able to derive analytical expressions for the basis functions and their corresponding eigenvalues. This will be much cheaper than computing the numerical approximation.
We want to only keep the eigenvalues (and corresponding eigenfunctions) of the Karhunen expansion that are above \epsilon = 10^{-16}. Everything else is obscured by numerical noise in the algorithms we feed the random field into.
Depending on the coefficient in the exponent of the function rho (in my MWE it was -5, but it has been as high as -200, depending on the simulation), it’s technically nonzero everywhere but that does not really matter when we can only resolve \epsilon. That was why we were making the matrix sparse for an approximation.
If for smaller coefficients where a sparse matrix approximation is a bad idea, what would be the best way to get around building the dense matrix if sparsity is no longer a good approximation?
If you want to set values below a certain threshold to zero, you must dismiss them before setting up the matrix, as otherwise the memory will be allocated, as is the case in the example above. The check if the value is below the threshold would have to be outside of rho, maybe like this:
from dolfin import *
from mshr import *
import scipy.sparse as sp
import numpy as np
from petsc4py import PETSc
mesh = generate_mesh(Circle(Point(0, 0), 1) - Circle(Point(0, 0), 0.8), 20)
V = FunctionSpace(mesh, "DG", 0)
def rho(x, y):
val = np.exp(-5 * (np.power(x[0] - y[0], 2) + np.power(x[1] - y[1], 2)) )
return val
xs = V.tabulate_dof_coordinates()
data = []
row_indices = []
col_indices = []
nnz = 0
for i in range(len(xs)):
for j in range(len(xs)):
val = rho(xs[i, :], xs[j, :])
if val > 1e-3:
nnz += 1
row_indices.append(i)
col_indices.append(j)
data.append(val)
C_intermediate = sp.csr_matrix( (data, (row_indices, col_indices)), shape=(V.dim(), V.dim()) )
C = PETScMatrix( PETSc.Mat().createAIJ(size = C_intermediate.shape, csr = (C_intermediate.indptr, C_intermediate.indices, C_intermediate.data)) )
print(nnz/V.dim()**2, nnz, C.nnz())