quickly parameter-dependent assembly

[AnyWayOne]

I wonder is it possible to implement quickly parameter-dependent assembly, e.g. the problem $-\triangledown(k(p)\triangledown(u))=1$, where $k(p)$ is a function related to the parameter $p$.
A normal implementation is:
##############
mesh = Mesh...
Vh = Basis...
@BilinearForm
def a(u,v,w);
return dot(w['kp']*grad(u), grad(v))
for k_p_i in k_p:
A_i = a.assemble(Vh, kp = Vh.interpolate(k_p_i))
####################
$A_i \in R^{n \times n}$
By using a loop, it executes n_p (the number of the samplr of the parameter $p$) assemblys, but it's too time consuming for a large $n_p$.
Is it possible to implement above loop efficiently at once and to obtain a tensor $A \in R^{n \times n \times n_p}$?
Thanks for any answer.

[kinnala]

There is no command to do them all at once. In principle it would be possible to have a "TrilinearForm" to assemble 3-tensors and then do multiple contractions to 2-tensors. However, I think the limitation is that there is no efficient data structure (similar to scipy.sparse) available for such sparse 3-tensors that would be suitable for FE assembly (sums multiple values during init).

One thing you could try is using Dask (and dask.bag) to parallelize the for-loop. Those computations do not depend on one another so if you have n_p computers, it should be pretty quick.

[AnyWayOne]

I strongly suggest adding the code into skfem, because it will be very efficient and convenient for parameter analysis, sensitivity analysis and uncertainty quantification, etc.

[AnyWayOne]

Currently, dense matrices are acceptable to my problems. It will become sparse again after a simple computation. It would be best if a sparse form can be directly obtained.

[kinnala]

Adding to the above code here is a snippet that I used:

def asm():
    m = MeshTri().refined(7)
    basis = Basis(m, ElementTriP1())
    from skfem.helpers import dot, grad
    out = TrilinearForm(lambda u, v, w, p: dot(w * grad(u), grad(v))).assemble(basis)
    import sparse
    arr = sparse.COO(np.array([out[1].flatten(), out[2].flatten(), out[3].flatten()]), out[0].flatten(), out[4], has_duplicates=True)
    k_p = np.random.rand(100, basis.N)
    res = sparse.tensordot(arr, k_p.T, axes=1, return_type=sparse.COO)

    A = []
    for mat in res.T:
        A.append(mat.to_scipy_sparse().tocsr())
    return A

It seems to run pretty quickly. However, I don't want to introduce an additional hard dependency to sparse because it depends on numba which is a rather heavy dependency and we try to keep it minimal. Therefore I think the best approach is to return something which is directly compatible with the constructor of sparse.COO.

[AnyWayOne]

It seems to be quick. I'll test it later on larger problems.
For the output of the sparse 3-tensor, generally there is no guarantee that the matrix with the third dimension is sparse, I think a possible format is to slice in the third dimension, i.e., A(:,:,i) is constructed as a sparse matrix, but A(i,:,:) or A(:,i,:) may be not a sparse matrix. This is just an extension of the concept of a two-dimensional sparse matrix. I have no idea about the general purpose of 3-tensors, this is only from the point of my problems.

[kinnala]

Continuing to track this is #733.