# Composite FE spaces + Parallel assembly with PETSc

**URL:** https://community.freefem.org/t/composite-fe-spaces-parallel-assembly-with-petsc/3336
**Category:** General Discussion
**Created:** [June 21, 2024, 4:44pm UTC](https://community.freefem.org/t/composite-fe-spaces-parallel-assembly-with-petsc/3336 "2024-06-21T16:44:01Z")
**Posts on this page:** 7
**Page:** 1

<div class="post-metadata">

### Author: ![pierreT](https://avatars.discourse-cdn.com/v4/letter/p/d07c76/32.png) [@pierreT](https://community.freefem.org/u/pierreT)
#### Post date: [June 21, 2024, 4:44pm UTC](https://community.freefem.org/t/composite-fe-spaces-parallel-assembly-with-petsc/3336/1 "2024-06-21T16:44:01Z")

</div>

Hello,  
I am a user of version 4.14, and I’ve run the test cases presented here:  
[https://doc.freefem.org/documentation/composite.html#with-petsc](https://doc.freefem.org/documentation/composite.html#with-petsc)  
I have a problem with parallel assembly + PETSc. Compared with the parallel assembly with MUMPS, I get a ratio of around 60 on matrix assembly time in favor of MUMPS over PETSc.  
I attach my scripts below.  
The way I run the scripts are “time ~/path\_to\_FreeFem/bin/ff-mpirun -np 4 ./test\_mumps.edp -wg”  
Calculating on a laptop running ubuntu 20.04.  
Thank you for help.  
Pierre

* * *

load “MUMPS”  
load “bem”

int nn = 200; // number of edges in each direction

mesh ThP = square(nn,nn,[2_pi_x,2_pi_y],flags=3); // Pressure mesh  
mesh ThU = trunc(ThP,1,split=2); // Velocity mesh

fespace Uh(ThU,[P1,P1]); // Velocity space  
fespace Ph(ThP,P1); // Pressure space

macro grad(u) [dx(u),dy(u)] //  
macro Grad(u1,u2) [grad(u1), grad(u2)] //  
macro div(u1,u2) (dx(u1)+dy(u2)) //

// definition of the boundary condition  
func g1 = sin(x)\*cos(y);  
func g2 = -cos(x)\*sin(y);

// definition of the right-hand side  
func f1 = 0;  
func f2 = -4\*cos(x)\*sin(y);

Uh [u1,u2],[v1,v2];  
Ph p,q;

fespace Xh=Uh\*Ph;

varf Stokes (\<[u1,u2],[p]\>, \<[v1,v2],[q]\>) = int2d(ThU)((Grad(u1,u2):Grad(v1,v2)))

- int2d(ThU)(-div(u1,u2)\*q -div(v1,v2)\*p)
- int2d(ThP)(-1e-10_p_q)
- int2d(ThU)([f1,f2]'\*[v1,v2])
- on(1,2,3,4, u1=g1, u2=g2);

real t1 = clock();  
matrix M = Stokes(Xh,Xh,solver=sparsesolver, master=-1);  
real t2 = clock();  
cout \<\< “init matrix time=” \<\< t2-t1 \<\< endl;

real[int] b = Stokes(0,Xh);  
real t3 = clock();  
cout \<\< “init rhs time=” \<\< t3-t2 \<\< endl;

real[int] sol = M^-1\*b;  
real t4 = clock();  
cout \<\< “solver time=” \<\< t4-t3 \<\< endl;

* * *

load “PETSc”  
load “bem”

int nn = 200; // number of edges in each direction

mesh ThP = square(nn,nn,[2_pi_x,2_pi_y],flags=3); // Pressure mesh  
mesh ThU = trunc(ThP,1,split=2); // Velocity mesh

fespace Uh(ThU,[P1,P1]); // Velocity space  
fespace Ph(ThP,P1); // Pressure space

macro grad(u) [dx(u),dy(u)] //  
macro Grad(u1,u2) [grad(u1), grad(u2)] //  
macro div(u1,u2) (dx(u1)+dy(u2)) //

// definition of the boundary condition  
func g1 = sin(x)\*cos(y);  
func g2 = -cos(x)\*sin(y);

// definition of the right-hand side  
func f1 = 0;  
func f2 = -4\*cos(x)\*sin(y);

Uh [u1,u2],[v1,v2];  
Ph p,q;

fespace Xh=Uh\*Ph;

varf Stokes (\<[u1,u2],[p]\>, \<[v1,v2],[q]\>) = int2d(ThU)((Grad(u1,u2):Grad(v1,v2)))

- int2d(ThU)(-div(u1,u2)\*q -div(v1,v2)\*p)
- int2d(ThP)(-1e-10_p_q)
- int2d(ThU)([f1,f2]'\*[v1,v2])
- on(1,2,3,4, u1=g1, u2=g2);

real t1 = clock();  
Mat M = Stokes(Xh,Xh);  
set(M, sparams = “-pc\_type lu”);  
real t2 = clock();  
cout \<\< “init matrix time=” \<\< t2-t1 \<\< endl;

real[int] b = Stokes(0,Xh);  
real t3 = clock();  
cout \<\< “init rhs time=” \<\< t3-t2 \<\< endl;

real[int] sol = M^-1\*b;  
real t4 = clock();  
cout \<\< “solver time=” \<\< t4-t3 \<\< endl;

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [June 21, 2024, 4:48pm UTC](https://community.freefem.org/t/composite-fe-spaces-parallel-assembly-with-petsc/3336/2 "2024-06-21T16:48:58Z")

</div>

Your code is not runnable, because it’s not in a code block, so we can’t help you. Furthermore, you don’t need composite fespaces for this particular problem.

---

<div class="post-metadata">

### Author: ![pierreT](https://avatars.discourse-cdn.com/v4/letter/p/d07c76/32.png) [@pierreT](https://community.freefem.org/u/pierreT)
#### Post date: [June 21, 2024, 4:57pm UTC](https://community.freefem.org/t/composite-fe-spaces-parallel-assembly-with-petsc/3336/3 "2024-06-21T16:57:06Z")

</div>

Sorry, here are code blocks. With MUMPS :

```auto
load "MUMPS"
load "bem"

int nn = 200; // number of edges in each direction

mesh ThP = square(nn,nn,[2*pi*x,2*pi*y],flags=3); // Pressure mesh
mesh ThU = trunc(ThP,1,split=2); // Velocity mesh

fespace Uh(ThU,[P1,P1]); // Velocity space
fespace Ph(ThP,P1); // Pressure space

macro grad(u) [dx(u),dy(u)] //
macro Grad(u1,u2) [grad(u1), grad(u2)] //
macro div(u1,u2) (dx(u1)+dy(u2)) //

// definition of the boundary condition
func g1 = sin(x)*cos(y);
func g2 = -cos(x)*sin(y);

// definition of the right-hand side
func f1 = 0;
func f2 = -4*cos(x)*sin(y);

Uh [u1,u2],[v1,v2];
Ph p,q;

fespace Xh=Uh*Ph;

varf Stokes (<[u1,u2],[p]>, <[v1,v2],[q]>) = int2d(ThU)((Grad(u1,u2):Grad(v1,v2)))
+ int2d(ThU)(-div(u1,u2)*q -div(v1,v2)*p)
+ int2d(ThP)(-1e-10*p*q)
+ int2d(ThU)([f1,f2]'*[v1,v2])
+ on(1,2,3,4, u1=g1, u2=g2);

real t1 = clock();
matrix M = Stokes(Xh,Xh,solver=sparsesolver, master=-1);
real t2 = clock();
cout << "init matrix time=" << t2-t1 << endl;

real[int] b = Stokes(0,Xh);
real t3 = clock();
cout << "init rhs time=" << t3-t2 << endl;

real[int] sol = M^-1*b;
real t4 = clock();
cout << "solver time=" << t4-t3 << endl;

```

and with PETSc :

```auto
load "PETSc"
load "bem"

int nn = 200; // number of edges in each direction

mesh ThP = square(nn,nn,[2*pi*x,2*pi*y],flags=3); // Pressure mesh
mesh ThU = trunc(ThP,1,split=2); // Velocity mesh

fespace Uh(ThU,[P1,P1]); // Velocity space
fespace Ph(ThP,P1); // Pressure space

macro grad(u) [dx(u),dy(u)] //
macro Grad(u1,u2) [grad(u1), grad(u2)] //
macro div(u1,u2) (dx(u1)+dy(u2)) //

// definition of the boundary condition
func g1 = sin(x)*cos(y);
func g2 = -cos(x)*sin(y);

// definition of the right-hand side
func f1 = 0;
func f2 = -4*cos(x)*sin(y);

Uh [u1,u2],[v1,v2];
Ph p,q;

fespace Xh=Uh*Ph;

varf Stokes (<[u1,u2],[p]>, <[v1,v2],[q]>) = int2d(ThU)((Grad(u1,u2):Grad(v1,v2)))
+ int2d(ThU)(-div(u1,u2)*q -div(v1,v2)*p)
+ int2d(ThP)(-1e-10*p*q)
+ int2d(ThU)([f1,f2]'*[v1,v2])
+ on(1,2,3,4, u1=g1, u2=g2);

real t1 = clock();
Mat M = Stokes(Xh,Xh);
set(M, sparams = "-pc_type lu");
real t2 = clock();
cout << "init matrix time=" << t2-t1 << endl;

real[int] b = Stokes(0,Xh);
real t3 = clock();
cout << "init rhs time=" << t3-t2 << endl;

real[int] sol = M^-1*b;
real t4 = clock();
cout << "solver time=" << t4-t3 << endl;

```

Thanks

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [June 21, 2024, 5:06pm UTC](https://community.freefem.org/t/composite-fe-spaces-parallel-assembly-with-petsc/3336/4 "2024-06-21T17:06:34Z")

</div>

The problem is in the line `Mat M = Stokes(Xh,Xh)`. I didn’t wrote that part of the code, so I can’t tell you why, but it’s highly inefficient, it’s not related to MUMPS. Again, you don’t really need composite fespaces for this, assembling `Mat M = Stokes(Xh,Xh)` with non-composite fespaces is just as fast as a standard assembly.

---

<div class="post-metadata">

### Author: ![phtournier](https://avatars.discourse-cdn.com/v4/letter/p/3da27b/32.png) [@phtournier](https://community.freefem.org/u/phtournier)
#### Post date: [June 22, 2024, 8:25pm UTC](https://community.freefem.org/t/composite-fe-spaces-parallel-assembly-with-petsc/3336/5 "2024-06-22T20:25:15Z")

</div>

Indeed, this is not intended behavior. Thank you for reporting the problem, I will look into it.

---

<div class="post-metadata">

### Author: ![phtournier](https://avatars.discourse-cdn.com/v4/letter/p/3da27b/32.png) [@phtournier](https://community.freefem.org/u/phtournier)
#### Post date: [June 22, 2024, 9:38pm UTC](https://community.freefem.org/t/composite-fe-spaces-parallel-assembly-with-petsc/3336/6 "2024-06-22T21:38:58Z")

</div>

This is fixed @ [fix computation of mapping for PETSc Mats with composite FEspaces · FreeFem/FreeFem-sources@4871f25 · GitHub](https://github.com/FreeFem/FreeFem-sources/commit/4871f258eb0e08c566a3cf0331404c5e377269ed)  
It will make it into the next release.

---

<div class="post-metadata">

### Author: ![pierreT](https://avatars.discourse-cdn.com/v4/letter/p/d07c76/32.png) [@pierreT](https://community.freefem.org/u/pierreT)
#### Post date: [June 23, 2024, 8:30pm UTC](https://community.freefem.org/t/composite-fe-spaces-parallel-assembly-with-petsc/3336/7 "2024-06-23T20:30:44Z")

</div>

@prj and @phtournier ,  
Thank you very much for your feedback. I’ll be testing the patch quickly.  
Pierre