Linear Equation from FFC generated from variational formulation

variational formulation
1

Hi FEniCS family,

I am working on on the following variational statement taken from
https://comet-fenics.readthedocs.io/en/latest/demo/periodic_homog_elas/periodic_homog_elas.html

Ve = VectorElement("CG", mesh.ufl_cell(), 2)
Re = VectorElement("R", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, MixedElement([Ve, Re]), constrained_domain=PeriodicBoundary(vertices, tolerance=1e-10))
V = FunctionSpace(mesh, Ve)

v_,lamb_ = TestFunctions(W)
dv, dlamb = TrialFunctions(W)
w = Function(W)
dx = Measure('dx')(subdomain_data=subdomains)

Eps = Constant(((0, 0), (0, 0)))
F = sum([inner(sigma(dv, i, Eps), eps(v_))*dx(i) for i in range(nphases)])
a, L = lhs(F), rhs(F)
a += dot(lamb_,dv)*dx + dot(dlamb,v_)*dx

For the variational statement a, L = lhs(F), rhs(F), the FEniCS find a linear equation (using FFC ) and after solve (a==L,w,[]) give output of fluctuation function, w for each node.

I wanted to find the linear system of equation generated by FEniCS (like Ax=b, with x=w being the fluctuation function). I want to use coefficient matrix A for my post processing.

I know that, FFC forms a C++ code from symbolic variational statement, and after solution gives back, function w in python format. So, there must be some linear equation which is generated to solve the variational form (a==L). Kindly tell me, how I can get access to the linear system of equation (coeff matrix, A being assembled over all elements of mesh, through guassian quadrature) and interpret it for usage.

Any help is greatly appreciated.

2

Try

A = assemble(a)
1 Like
3

Thanks @dokken for your response.
I am getting as object type A, and I am not able to use or extract any information from A.

A=assemble(a)
print(A)

<dolfin.cpp.la.Matrix object at 0x7ff36b8cc950>

Also, for general view, when we got this hexadecimal location. How can we get the object data/ value?

4

Following Transform dolfinx assemble matrix to numpy array - #11 by Samuel_Groth, but adapted to legacy DOLFIN, I would do something along the lines of

from scipy.sparse import csr_matrix
ai, aj, av = as_backend_type(A).mat().getValuesCSR()
Asp = csr_matrix((av, aj, ai))

to get a SCIPY CSR matrix.

1 Like
5

Thanks @dokken for your response.

A=assemble(a)
from scipy.sparse import csr_matrix
ai, aj, av = as_backend_type(A).mat().getValuesCSR()
Asp = csr_matrix((av, aj, ai))
Asp
<19851x19851 sparse matrix of type '<class 'numpy.float64'>
	with 535905 stored elements in Compressed Sparse Row format>

Can you tell me little more about ai,aj,av. I read about csr_matrix(data,(row,column)) and I was trying to relate with ai, aj and av based on their size.

len(ai)
19852
len(aj)
535905
len(av)
535905

Asp = csr_matrix((av, aj, ai)).toarray()
# giving Asp as matrix of size 19851 X 19851

But, I am getting difficulty of what they try to interpret. I tried with mesh data.

mesh.num_vertices()
6771
mesh.num_cells()
13232

Can you tell me how Asp is interpreted with respect to form a (bilinear variation form used in Assemble(a))?
Also, in as_backend_type, Does Asp gives coefficent matrix (K) from Kx=f (say in general with x as fluctuation function to be evaluated) ?
and can we evaluate the vector f, since, as_backend_type also return a vector.

Thanks for your help.

6

Please read the scipy docs: scipy.sparse.csr_matrix — SciPy v1.10.1 Manual

You are using a second order Lagrange space, with the addition of a real space dof. You cannot relate those to mesh entities.

b =assemble(L).get_local()
x = w.vector().get_local()

aligns as Ax=b.

Note that his would only work in serial, as numpy/scipy based matrices do not have parallel support (as in being distributed).

1 Like
7

Thanks @dokken for your response. It helped me to identify the linear system.

The Asp (coefficient matrix) formed is of very large size even with linear elements.

A=assemble(a)
from scipy.sparse import csr_matrix
ai, aj, av = as_backend_type(A).mat().getValuesCSR()
Asp = csr_matrix((av, aj, ai)).toarray()

# Asp: 19851x19851 sparse matrix of type '<class 'numpy.float64'>
	#with 535905 stored elements in Compressed Sparse Row format>

Is there any suitable method/ function in FEniCS to invert Asp (coefficient) matrix efficiently and accurately ? (Mainly the numerical method)

I tried the following:

from scipy.sparse import csc_matrix
from scipy.sparse.linalg import inv
A = csc_matrix(Asp)
inv(A)

But, it crashes the kernel

8

The whole point of using dolfin is that it invert this matrix (using petsc), and gives you a solution to the problem. You will never want to store the actual inverse, as this is a dense system

if you use a direct solver (LUSolver), one can store the LU factorization, which makes the solve step quick. See for instance: https://fenicsproject.org/qa/6045/repeated-use-of-matrix-with-lu-decomposition/

https://fenicsproject.org/pub/tutorial/html/._ftut1017.html

1 Like
9

Thanks @dokken for your response. I have few doubts:

  1. The linear system which FFC gives is of form Aw=b, where FFC performs some calculation to get A^{-1} to solve for
    w=A^{-1}b.
    Is there any way to work on FFC code / PETSc solver to get that A^{-1} ?

  2. Can I know the intermediate linear form which comes just after minimization of variational statement (a), which later FFC convert to Aw=b form.

For example : \int \:\sigma(w)\:\ grad(v)\:dvol \: where
w : trial function (fluctuation function), v: test function, dvol: volume

gives the minimization output of U\:w\:=\:P\:\epsilon ....(1) analytically, where epsilon is given macro-strain [6X1] vector given as Eps.
[in periodic homogenization code (FEniCs example)].

However, the output given by FFC through sparse csr is
A w =b;
instead of initial U\:w\:=\:P\:\epsilon , where
U= A and P \epsilon=b

How can I get the intermediate term P ?
Can we know the form generated by FFC eqn …(1) instead of direct A w =b ?

Any help is appreciated.

10

Lets make some things clear here:

  1. FFC does not create the in inverse. FFC takes a ufl form, and generates code for integrating over a cell/facet to return the local tensor (scalar, vector, matrix) for a single cell/facet.
  2. A= assemble(a) generates the sparse matrix by calling the kernels from step 1 for all cells/facets specified in the form.
  3. this matrix is sent to a linear algebra backend, such as PETSc, that uses a direct or iterative solver to get the solution to the problem Ax=b. At no point is the inverse matrix actually formed. See: KSP: Linear System Solvers — PETSc 3.18.5 documentation
    For instance it states

To solve a linear system with a direct solver (currently supported by PETSc for sequential matrices, and by several external solvers through PETSc interfaces, see Using External Linear Solvers) one may use the options -ksp_type preonly (or the equivalent -ksp_type none) -pc_type lu (see below).

For instance you can use the LU factorizer UMFPACK (MATSOLVERUMFPACK — PETSc 3.18.5 documentation) or mumps (MATSOLVERMUMPS — PETSc 3.18.5 documentation)

An LU factorization never forms the full inverse, see: LU decomposition - Wikipedia

1 Like
11

Thanks @dokken for your response.

I am facing problem with integrating float type tensor over all the cells.

Em = 50e3
num = 0.2
Er = 210e3
nur = 0.3
material_parameters = [(Em, num), (Er, nur)]
nphases = len(material_parameters)

def Com(i):   
    E, nu = material_parameters[i]
    lmbda = E*nu/(1+nu)/(1-2*nu)
    mu = E/2/(1+nu)
    C1=lmbda+2*mu
    return Constant([(C1,lmbda,lmbda,0,0,0),(lmbda,C1,lmbda,0,0,0),(lmbda,lmbda,C1,0,0,0),(0,0,0,mu,0,0),(0,0,0,0,mu,0),(0,0,0,0,0,mu)])
   

D= np.zeros((6,6))
for j in range(6):
    for k in range(6):
       D[j,k]= assemble(sum([(Com(i)[j,k])*dx(i) for i in range(nphases)]))/vol

This integral is missing an integration domain.
---------------------------------------------------------------------------
UFLException                              Traceback (most recent call last)
/tmp/ipykernel_455/2734975225.py in <module>
      2 for j in range(6):
      3     for k in range(6):
----> 4        D[j,k]= assemble(sum([(Com(i)[j,k])*dx(i) for i in range(nphases)]))/vol

/tmp/ipykernel_455/2734975225.py in <listcomp>(.0)
      2 for j in range(6):
      3     for k in range(6):
----> 4        D[j,k]= assemble(sum([(Com(i)[j,k])*dx(i) for i in range(nphases)]))/vol

/usr/lib/python3/dist-packages/ufl/measure.py in __rmul__(self, integrand)
    430                 domain, = domains
    431             elif len(domains) == 0:
--> 432                 error("This integral is missing an integration domain.")
    433             else:
    434                 error("Multiple domains found, making the choice of integration domain ambiguous.")

/usr/lib/python3/dist-packages/ufl/log.py in error(self, *message)
    156         "Write error message and raise an exception."
    157         self._log.error(*message)
--> 158         raise self._exception_type(self._format_raw(*message))
    159 
    160     def begin(self, *message):

UFLException: This integral is missing an integration domain.

Can you help me to compute the assemble command ?

I am getting error with multiple attempts like changing Com function
return as_tensor([(C1,lmbda,lmbda,0,0,0),(lmbda,C1,lmbda,0,0,0),(lmbda,lmbda,C1,0,0,0),(0,0,0,mu,0,0),(0,0,0,0,mu,0),(0,0,0,0,0,mu)])
or
return ([(C1,lmbda,lmbda,0,0,0),(lmbda,C1,lmbda,0,0,0),(lmbda,lmbda,C1,0,0,0),(0,0,0,mu,0,0),(0,0,0,0,mu,0),(0,0,0,0,0,mu)])

Any help is greatly appreciated.

12

The code you have posted is not complete, as you have not defined dx and it is unclear what dx(i) is supposed to represent.
Please see:Read before posting: How do I get my question answered? - #4 by nate for examples following the forum guidelines

13

dx is taken from the above step and subdomain is similar to the periodic homogenization example.

from __future__ import print_function
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
%matplotlib notebook

a = 1.         # unit cell width
b = sqrt(3.)/2. # unit cell height
c = 0.5        # horizontal offset of top boundary
R = 0.2        # inclusion radius
vol = a*b      # unit cell volume
# we define the unit cell vertices coordinates for later use
vertices = np.array([[0, 0.],
                     [a, 0.],
                     [a+c, b],
                     [c, b]])
fname = "hexag_incl"
mesh = Mesh(fname + ".xml")
subdomains = MeshFunction("size_t", mesh, fname + "_physical_region.xml")
facets = MeshFunction("size_t", mesh, fname + "_facet_region.xml")
plt.figure()
plot(subdomains)
plt.show()

I apologize for the incomplete step.

14

Then change

To
dx = Measure('dx')(domain=mesh, subdomain_data=subdomains)

1 Like