Partitioning the stiffness matrix without partitioning the mesh for PETSc

General Discussion
1

Hi FreeFEM developers!

I thank you for having integrated PETSc with FreeFEM!
This has resulted in an easy implementation of large scale topology optimization.

I am attempting to implement a contact mechanics solver in parallel (for topology optimization).
Implemented via fixed point and penalization, the contact problem reads
(K + M(U_n))U_{n+1} = F
where U_n is the solution iterate, K is the linear elasticity matrix, F is the external force and M(U_n) is the contact penalization term.
To simplify the construction of M(U_n), I work with co-incident contact surfaces (in 2D and in 3D).
In this case, it suffices to determine indices of the coincident vertices, and the corresponding pairing matrix.
Using the pairing matrix, I construct the matrix M(U_n) by specifying values at certain indices (ones corresponding to the contact surface).

For ​PETSc, I need to create subdomains with ghost elements. This complicates the creation of the pairing matrix and the construction of the matrix M(U_n).
Is there a possibility to simply create the global matrix M(U_n) in sequential and then feed it to a macro like createMat?
If not, is there any other solution?

2

For ​PETSc, I need to create subdomains with ghost elements. This complicates the creation of the pairing matrix and the construction of the matrix M(U_n).

How? Do you have a concrete example to share?

Is there a possibility to simply create the global matrix M(U_n) in sequential

Well, if you know how to compute M(U_n) in sequential, then you could just restrict the global matrix on each subdomain. But I would highly advise against this solution.

3

In order to explain, I have attached an image.
subdomain_crack.pdf (7.0 KB)
In it, we see the crack is represented by \Gamma^+, \Gamma^- and there are three subdomains.
For each vertex on \Gamma^+, I need to find its nearest neighbor vertex on \Gamma^-.
In absence of subdomain, this is really easy (I can attach a code, if it helps to understand better).
In the presence of subdomains, I have to create a loop, checking for the neighbor vertex in every subdomain. Which is doable, but a bit messy.

The matrix M(U_N) has the structure as shown in the pdf below.
matrix_structure.pdf (7.4 KB)
In it, there the sub-matrix C_11 is easy to construct (simply using varf). But, C_12 ought to be constructed manually.
Each entry of matrix C_12, say v_ij is constructed by finding the value of v_kj (where index k is the nearest neighbor of index i). Doing this ensuing that the partition of unity is not violated (as I shall be looping on every subdomain containing \Gamma^+) seems to be a lot of work.

Given that I have already coded all of the above operations in sequential, I thought that it would be much easier if I could avoid re-coding everything for subdomains.

I hope I could explain my problem well. If needed I can provide more details.

Your suggestion is to create the interpolation matrix between the global and the local FE spaces, and then to multiply the matrix M(U_N) with the interpolation matrix right? Or even simpler?
Nice idea. Did not think of that. But why do you not recommend it?

4

Your suggestion is to create the interpolation matrix between the global and the local FE spaces

I never said to create the interpolation matrix, which will be painfully slow. You can restrict a global matrix using N2O + restrict.
Something like:

matrix LocalC12 = GlobalC12(LocalRows, LocalCols);
Mat GlobalC12(C11, C22, LocalC12); // rectangular distributed Mat

But why do you not recommend it?

It is sequential by nature, as soon as you start looking at large 3D meshes, it will probably become the bottleneck of your method. Maybe not, I’d need to see your implementation first.

5

Okay, my bad.
I understand that even if its a sparse matrix multiplication, it will be slow.

I did not completely understand this solution.
First I create the global matrix C12 and distribute the mesh using buildDmesh and n2o.
And then I use the two commands you suggested, followed by createMat(Th, GlobalC12, Pk)?

Here is the implementation in sequential:
(after creating the pairing vector, denoted by Pairing)

6

First I create the global matrix C12 and distribute the mesh using buildDmesh and n2o.
And then I use the two commands you suggested, followed by createMat(Th, GlobalC12, Pk)?

You said that you have two meshes \Gamma^+ and \Gamma^-, didn’t you? So you need to call twice buildDmesh with each mesh, then createMat to build the diagonal blocks. Then, you’ll be able to create a rectangular block following the distribution of both diagonal blocks using the local restricted GlobalC12.

But, from looking at your script (which is indeed quite slow I guess since you are doing this in FreeFEM instead of in a small C++ plugin), it looks like you need the full C11 as well to build C12. So, if you really want to do this, I guess there is no other option than to build the global M(U_n) and then restrict using N2O. Since M(U_n) is square, I guess it is as easy as matrix LocMUn = GlobalMUn(LocalUnknowns, LocalUnknowns); Mat MUn(A, LocMUn); where A is a distributed Mat following the domain decomposition of \Gamma.

7

Here is a small example in case it is not clear for you.

include "macro_ddm.idp"
load "PETSc"

mesh Th = square(20, 20);
mesh ThGlobal = Th;

int[int] n2o;
macro ThN2O()n2o//

Mat A;
createMat(Th, A, P2)

varf vPb(u, v) = int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v)) + on(1,2,3, u = 1.0) + int2d(Th)(v);
fespace Vh(Th, P2);
A = vPb(Vh, Vh, tgv = -1);
set(A, sparams = "-pc_type lu");
real[int] rhs = vPb(0, Vh, tgv = -1);
Vh sol;
sol[] = A^-1 * rhs;
plotD(Th, sol, cmm = "Solution from distributed assembly")
fespace VhGlobal(ThGlobal, P2);
int[int] localUnknowns = restrict(Vh, VhGlobal, n2o);
mesh ThBackup = Th;
Th = ThGlobal;
matrix global = vPb(VhGlobal, VhGlobal);
matrix local = global(localUnknowns, localUnknowns);
Th = ThBackup;
A = local;
sol[] = A^-1 * rhs;
plotD(Th, sol, cmm = "Solution from centralized assembly")
8

No, I have two interfaces representing the same surface. The mesh is the same.
On FreeFEM, given that one cannot create coincident surfaces in the same mesh, I create a mesh with some gap and then use movemesh to close the gap to zero.
So I guess I need to call buildDmesh only once.

Yes, a C++ plugin would have been much faster!

9

So I guess I need to call buildDmesh only once.

OK, so I guess you can just copy/paste what I wrote just above, and you should be good to go.
As long as you backup the global mesh and know how to go from local to global and vice versa, anything is possible :smile:

10

With this script, it is very clear.
Thank you so much!
This solves my problem!! :man_dancing: :man_dancing: :man_dancing:

11

I totally agree with you!
The restriction operator with the backup is very very useful!