# Exporting Matrices of AS method for Poisson Equation

**URL:** <https://community.freefem.org/t/exporting-matrices-of-as-method-for-poisson-equation/526>\
**Category:** General Discussion\
**Created:** [July 30, 2020, 8:02pm UTC](https://community.freefem.org/t/exporting-matrices-of-as-method-for-poisson-equation/526 "2020-07-30T20:02:39Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![gernot](https://avatars.discourse-cdn.com/v4/letter/g/dfb087/32.png) [@gernot](https://community.freefem.org/u/gernot)\
**Post date:** [July 30, 2020, 8:02pm UTC](https://community.freefem.org/t/exporting-matrices-of-as-method-for-poisson-equation/526/1 "2020-07-30T20:02:39Z")

</div>

I want to decompose a 3D mesh (sphere) into multiple subdomains and export the stiffness matrices / load vectors to solve the problem with an external CG solver.  
Unfortunately it seems that the local matrices are non-symmetric.

```
// ff-mpirun -np 4 sphere2.edp -wg -ffddm_overlap 2 -ffddm_schwarz_method asm

macro dimension 3 // EOM

verbosity = 0;

load "msh3"
load "tetgen"
// load "medit"
// load "iovtk"

include "MeshSurface.idp"
include "ffddm.idp"

real radius = 5;
real hh = radius * 0.1;
meshS ThS = Sphere(radius, hh, 1, 0., 0., 0., 1);
// meshS ThS = Sphere(radius, 500., 500., 500., hh, 1, 1); // laut Doku
// plot(ThS,wait=1);

real[int] domain =[0.,0.,0.,1,hh^3]; // x, y, z, attrib, max_vol
mesh3 Th3 = tetg(ThS, switch="pqaAAQYY", nbofregions=1, regionlist=domain);
// savevtk("sphere.vtu", Th3);

// Mesh refinement
cout << "Num vertices coarse: " << Th3.nv << ", hmax = " << Th3.hmax << endl;
// medit("Coarse", Th3);

int refineSteps = 0;
for (int i = 0; i < refineSteps; i++)
{
    Th3 = trunc(Th3, 1, split=2);
    cout << "Num vertices fine: " << Th3.nv << ", hmax = " << Th3.hmax << endl;
}
// medit("Fine", Th3);

func bc = x^3*z^2*sin(y);

macro Grad3(u) [dx(u),dy(u),dz(u)] // EOM

macro Varf(varfName, meshName, VhName)
    varf varfName(u,v) = int3d(meshName)(Grad3(u)'* Grad3(v)) + on(1, u=bc);
// EOM

// Domain decomposition
ffddmbuildDmesh( LapMesh , Th3 , mpiCommWorld )
// medit("sphere" + mpirank, LapMeshThi);
int[int] ll = labels(LapMeshThi);
meshS surf = extract(LapMeshThi, label=ll);
cout << "Rank " << mpirank << ": Num vertices " << LapMeshThi.nv << ", extracted vertices: " << surf.nv << endl;
// savevtk("Th" + mpirank + ".vtk", LapMeshThi);

macro def(u) u// EOM
macro init(u) u// EOM
ffddmbuildDfespace(LapFE, LapMesh, real, def, init, P1)
ffddmsetupOperator( PB , LapFE , Varf )

func f = -6*x*z^2*sin(y)+x^3*z^2*sin(y) -2*x^3*sin(y);
macro myVarfrhs(varfName, meshName, VhName)
    varf varfName(u,v) = int3d(meshName)(f*v) + on(1, u=bc); // f + Laplace g
// EOM

real[int] rhsi(0);
ffddmbuildrhs(PB, myVarfrhs, rhsi)

if (mpisize == 1)
{
    ofstream fout("A_global.txt");
    fout << PBaRd[0] << endl;
}
else
{
    ofstream fout("A_local" + mpirank + ".txt");
    fout << PBaRd[mpirank] << endl;
}

if (mpisize == 1)
{
    ofstream fout("R_global.txt");
    fout << LapFERih[0] << endl;
}
else
{
    ofstream fout("R_local" + mpirank + ".txt");
    fout << LapFERih[mpirank] << endl;
}

if (mpisize == 1)
{
    ofstream fout("rhs_global.txt");
    fout << rhsi << endl;
}
else
{
    ofstream fout("rhs_local" + mpirank + ".txt");
    fout << rhsi << endl;
}

if (mpisize == 1)
{
    ofstream fout("PoU_global.txt");
    fout << LapFEDih[0] << endl;
}
else
{
    ofstream fout("PoU_local" + mpirank + ".txt");
    fout << LapFEDih[mpirank] << 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:** [July 30, 2020, 8:18pm UTC](https://community.freefem.org/t/exporting-matrices-of-as-method-for-poisson-equation/526/2 "2020-07-30T20:18:44Z")

</div>

You can do that with PETSc quite trivially, with `-ksp_type cg -pc_type asm`. The interface is much more concise as well, see [this](http://jolivet.perso.enseeiht.fr/FreeFem-tutorial/section_8/example3.edp.html) example.

---

<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:** [July 30, 2020, 8:31pm UTC](https://community.freefem.org/t/exporting-matrices-of-as-method-for-poisson-equation/526/3 "2020-07-30T20:31:00Z")

</div>

Here is the script using PETSc CG + ASM. You can also export the matrix and precondition it externally using the additional command line argument ` -mat_view ascii:Poisson`.

```auto
macro dimension 3 // EOM

verbosity = 0;

load "msh3"
load "tetgen"
// load "medit"
// load "iovtk"

include "MeshSurface.idp"
load "PETSc"
macro def(u) u// EOM
macro init(u) u// EOM
include "macro_ddm.idp"

real radius = 5;
real hh = radius * 0.1;
meshS ThS = Sphere(radius, hh, 1, 0., 0., 0., 1);
// meshS ThS = Sphere(radius, 500., 500., 500., hh, 1, 1); // laut Doku
// plot(ThS,wait=1);

real[int] domain =[0.,0.,0.,1,hh^3]; // x, y, z, attrib, max_vol
mesh3 Th3 = tetg(ThS, switch="pqaAAQYY", nbofregions=1, regionlist=domain);
// savevtk("sphere.vtu", Th3);

// Mesh refinement
cout << "Num vertices coarse: " << Th3.nv << ", hmax = " << Th3.hmax << endl;
// medit("Coarse", Th3);

int refineSteps = 0;
for (int i = 0; i < refineSteps; i++)
{
    Th3 = trunc(Th3, 1, split=2);
    cout << "Num vertices fine: " << Th3.nv << ", hmax = " << Th3.hmax << endl;
}
// medit("Fine", Th3);

func bc = x^3*z^2*sin(y);

macro Grad3(u) [dx(u),dy(u),dz(u)] // EOM

// Domain decomposition
buildDmesh( Th3 )
// medit("sphere" + mpirank, Th3);
int[int] ll = labels(Th3);
meshS surf = extract(Th3, label=ll);
cout << "Rank " << mpirank << ": Num vertices " << Th3.nv << ", extracted vertices: " << surf.nv << endl;
// savevtk("Th" + mpirank + ".vtk", Th3);

Mat A;
createMat( Th3, A, P1 )

func f = -6*x*z^2*sin(y)+x^3*z^2*sin(y) -2*x^3*sin(y);
    varf Poisson(u,v) = int3d(Th3)(Grad3(u)'* Grad3(v)) + int3d(Th3)(f*v) + on(1, u=bc); // f + Laplace g
fespace Vh(Th3, P1);
A = Poisson(Vh, Vh, sym = 1);
real[int] rhs = Poisson(0, Vh);
set(A, sparams = "-pc_type asm -pc_asm_overlap 2 -pc_asm_type basic -sub_pc_type cholesky -ksp_type cg -ksp_monitor -ksp_view");
real[int] sol = A^-1 * rhs;

```

---

<div class="post-metadata">

**Author:** ![gernot](https://avatars.discourse-cdn.com/v4/letter/g/dfb087/32.png) [@gernot](https://community.freefem.org/u/gernot)\
**Post date:** [July 31, 2020, 8:21am UTC](https://community.freefem.org/t/exporting-matrices-of-as-method-for-poisson-equation/526/4 "2020-07-31T08:21:18Z")

</div>

Thank you for your script.

I have run it with “ff-mpirun -n 4 sphere3.edp -mat\_view ascii:Poisson” and I get a “Poisson” textfile containing a matrix description.  
But it seems it is only one partition? At least it says “Mat Object: 1 MPI processes”.

How can I get all the submatrices? Is there a way to visualize the subdomains?

---

<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:** [July 31, 2020, 8:26am UTC](https://community.freefem.org/t/exporting-matrices-of-as-method-for-poisson-equation/526/5 "2020-07-31T08:26:03Z")

</div>

You are correct, it’s the global matrix. You are free to partition this matrix how you seem appropriate in your external solver if you don’t want to use PETSc CG.

---

<div class="post-metadata">

**Author:** ![gernot](https://avatars.discourse-cdn.com/v4/letter/g/dfb087/32.png) [@gernot](https://community.freefem.org/u/gernot)\
**Post date:** [July 31, 2020, 8:28am UTC](https://community.freefem.org/t/exporting-matrices-of-as-method-for-poisson-equation/526/6 "2020-07-31T08:28:51Z")

</div>

Is there a way to export the local matrices from the ASM as well?

I would like to implement the whole AS+CG algorithm and then later compare it to the solution of FreeFem. I just need some program to do the meshing/partitioning/assembly.

---

<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:** [July 31, 2020, 8:39am UTC](https://community.freefem.org/t/exporting-matrices-of-as-method-for-poisson-equation/526/7 "2020-07-31T08:39:33Z")

</div>

Here is what you can do:

1. replace `sym = 1` by `sym = 0`
2. replace the current `sparams` by:

```auto
set(A, sparams = "-pc_type gasm -pc_gasm_overlap 2 -pc_gasm_type basic -sub_pc_type cholesky -ksp_type cg -ksp_monitor -ksp_view -pc_gasm_view_subdomains -ksp_view_mat ascii:Poisson:ascii_matlab -ksp_view_rhs ascii:rhs:ascii_matlab -sub_ksp_view_mat ascii:Poisson_"+mpirank+":ascii_matlab");

```

This will export a `Poisson` file which is the global system, and `Poisson_rank` for each local subdomain matrices. In the KSP view, you’ll see each subdomain indices (inner means without overlap, outer means with overlap). It also saves the global right-hand side.

---

<div class="post-metadata">

**Author:** ![gernot](https://avatars.discourse-cdn.com/v4/letter/g/dfb087/32.png) [@gernot](https://community.freefem.org/u/gernot)\
**Post date:** [July 31, 2020, 12:23pm UTC](https://community.freefem.org/t/exporting-matrices-of-as-method-for-poisson-equation/526/8 "2020-07-31T12:23:34Z")

</div>

Thank you very much, that was really helpful.

Can I export the restriction matrices and the decomposition of the one as well?

---

<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:** [July 31, 2020, 5:32pm UTC](https://community.freefem.org/t/exporting-matrices-of-as-method-for-poisson-equation/526/9 "2020-07-31T17:32:22Z")

</div>

The restriction is the list of indices in the `Outer subdomain:` that you get with `-ksp_view`. You can define in your external solver a matrix of dimension (size of local matrix, that you get in the `-ksp_view` as well) x (size of global matrix, in your case 5654). Then, you just fill with ones for all lines at each column index given in the `-ksp_view`.