# Global solution with PETSc

**URL:** <https://community.freefem.org/t/global-solution-with-petsc/239>\
**Category:** General Discussion\
**Created:** [February 6, 2020, 1:23pm UTC](https://community.freefem.org/t/global-solution-with-petsc/239 "2020-02-06T13:23:18Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![henkel](https://avatars.discourse-cdn.com/v4/letter/h/b77776/32.png) [@henkel](https://community.freefem.org/u/henkel)\
**Post date:** [February 6, 2020, 1:23pm UTC](https://community.freefem.org/t/global-solution-with-petsc/239/1 "2020-02-06T13:23:18Z")

</div>

For verification purposes I would like to compare a solution computed in parallel with the sequentially computed one.

As a basis for an minimal working example I took an [example from Pierre Jolivet’s tutorial](http://jolivet.perso.enseeiht.fr/FreeFem-tutorial/section_8/example4.edp) and tried to reduce the distributed local solutions to a global one. I tried to accomplish this in the same way as it was done in the [maxwell-3d-PETSc.edp](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/maxwell-3d-PETSc.edp) example:

```
load "PETSc"
load "medit"

macro dimension()2// EOM
include "macro_ddm.idp"

real c = 6.25;
int N = 80;

mesh M;
{
    border a(t=0,1) { x=t; y=0; label=1; };
    border b(t=0,1) { x=1; y=t; label=1; };
    border c(t=0,1) { x=1-t; y=1; label=1; };
    border d(t=0,1) { x=0; y=1-t; label=1; };
    M = buildmesh(a(N)+b(N)+c(N)+d(N));
}

mesh ThPlt = M;

load "Element_P3"
func Pk = P3;
Mat J;

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

createMat(M, J, Pk)

fespace Vh(M, Pk);
Vh u;
func BC = cos(pi*x)*cos(pi*y);
varf vInit(w, v) = on(1, w = BC);
varf vJ(w, v) = int2d(M)(dx(w)*dx(v) + dy(w)*dy(v) - c*exp(u)*w*v) + on(1, w = 0);
varf vRes(w, v) = int2d(M)(dx(u)*dx(v) + dy(u)*dy(v) - c*exp(u)*v) + on(1, w = u);
int it = 1;
real tol = 1.0e-6;
u[] = vInit(0, Vh, tgv = -1);
real[int] b = u[];
for(int i = 1; i < it; ++i) {
    real[int] v = vRes(0, Vh, tgv = -1);
    v -= b;
    real res = sqrt(J(v, v));
    if(mpirank == 0)
        cout << "Norm of the residual = " << res << endl;
    if(res < tol)
        break;
    matrix Loc = vJ(Vh, Vh, tgv = -1);
    J = Loc;
    real[int] inc = J^-1 * v;
    u[] -= inc;
    plotD(M, u, cmm = "Solution at iteration #" + i);
}
plotD(M, u, cmm = "Final solution");

if(!NoGraphicWindow) {
    real[int] sol;
    changeNumbering(J, u[], sol);
    changeNumbering(J, u[], sol, inverse = true);
    fespace WhPlt(ThPlt, Pk);
    WhPlt plt;
    WhPlt reduce;
    cout << mpirank << " " << u[].n << " " << sol.n << " " << Vh.ndof << " " << WhPlt.ndof << " " <<endl;
    int[int] rest = restrict(Vh, WhPlt, n2o);
    for[i, v : rest] plt[][v] = sol[i];
    mpiReduce(plt[], reduce[], processor(0, mpiCommWorld), mpiSUM);
    if(mpirank == 0)
        medit("Global solution", ThPlt, reduce); // To keep this example short I kept medit here instead of saving the solution or calculating the difference to a reference solution
}

```

But the restrict(?) results in a segfault:  
`[1]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range`

I would appreciate it if someone would point out my mistake.

---

<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:** [February 6, 2020, 1:46pm UTC](https://community.freefem.org/t/global-solution-with-petsc/239/2 "2020-02-06T13:46:22Z")

</div>

You made two small mistakes which explain the segfault. Here is the fix:

1. `macro ThN2O()n2o// EOM` should read `macro MN2O()n2o// EOM` because your mesh name is `M`, not `Th` as in `maxwell-3d-PETSc.edp`
2. `changeNumbering(J, u[], sol);` goes from FreeFEM to PETSc numbering, then `changeNumbering(J, u[], sol, inverse = true);` goes from PETSc to FreeFEM numbering. You want to look at the values of `u[]`, not `sol`. `for[i, v : rest] plt[][v] = sol[i];` is wrong and should indeed be `for[i, v : rest] plt[][v] = u[][i];`

Let me know if I can be of any further assistance.

---

<div class="post-metadata">

**Author:** ![macro2micro](https://avatars.discourse-cdn.com/v4/letter/m/13edae/32.png) [@macro2micro](https://community.freefem.org/u/macro2micro)\
**Post date:** [April 9, 2020, 10:08am UTC](https://community.freefem.org/t/global-solution-with-petsc/239/3 "2020-04-09T10:08:46Z")

</div>

I can’t reply to you in that way. How can i touch with you?

---

<div class="post-metadata">

**Author:** ![macro2micro](https://avatars.discourse-cdn.com/v4/letter/m/13edae/32.png) [@macro2micro](https://community.freefem.org/u/macro2micro)\
**Post date:** [April 9, 2020, 10:09am UTC](https://community.freefem.org/t/global-solution-with-petsc/239/4 "2020-04-09T10:09:26Z")

</div>

![2222222](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/d3a6170fa4936d8ba8e9307925f4e916a658a8bb.png)

---

<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:** [April 9, 2020, 10:20am UTC](https://community.freefem.org/t/global-solution-with-petsc/239/5 "2020-04-09T10:20:29Z")

</div>

Please do not spam other threads…

---

<div class="post-metadata">

**Author:** ![jgonzalez49](https://avatars.discourse-cdn.com/v4/letter/j/f05b48/32.png) [@jgonzalez49](https://community.freefem.org/u/jgonzalez49)\
**Post date:** [August 6, 2020, 10:14pm UTC](https://community.freefem.org/t/global-solution-with-petsc/239/6 "2020-08-06T22:14:22Z")

</div>

I seem to be having the same issue with the segmentation fault. I have taken the [stokes 3d PETSc example](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/stokes-3d-PETSc.edp) and adapted it as I understood it from the original example above and the corrections from @prj. Here is the resulting edp file:

[stokes-3d-PETScagain.edp](https://community.freefem.org/uploads/short-url/jrKwrRaBMfhVLT05uKoEHyqPKiw.edp) (2.1 KB)

By systematically removing lines, I too seem to find the error only once the line

```
int[int] rest = restrict(Whbackup, Wh, n2o);

```

is included. As best I can tell, the two things included in the corrections are addressed since my original mesh name is `Th` for the macro line, and I use the values of `u[]` in the for loop. I appreciate any suggestions on fixing this issue.

---

<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:** [August 7, 2020, 5:46am UTC](https://community.freefem.org/t/global-solution-with-petsc/239/7 "2020-08-07T05:46:41Z")

</div>

You are using `buildMat`, you need to use `createMat` instead. With the following line instead, I don’t get any segmentation fault.

```auto
createMat(Th, A, Pk)

```

---

<div class="post-metadata">

**Author:** ![jgonzalez49](https://avatars.discourse-cdn.com/v4/letter/j/f05b48/32.png) [@jgonzalez49](https://community.freefem.org/u/jgonzalez49)\
**Post date:** [August 7, 2020, 2:25pm UTC](https://community.freefem.org/t/global-solution-with-petsc/239/8 "2020-08-07T14:25:28Z")

</div>

That works perfectly, thank you. I’ve been reviewing the other active thread [here](https://community.freefem.org/t/saving-and-loading-global-distibuted-arrays/493) and now see that the createMat command is common to the two examples you had mentioned at the top that specifically gather the results at the end of the program.

As a followup question (maybe better moved to the below mentioned thread if more pertinent), I notice that the file output of this now matches the [description of the output files of another user](https://community.freefem.org/t/plotmpi-how-to-get-plot-display-for-parallelization-hpddm/162/5) when they were using an older version of one of the examples and generated a single PVD file and four VTU files. The suggestion you gave made use of a [revised version of that example](https://github.com/FreeFem/FreeFem-sources/commit/ba61bcbfd00e10a44e750ba15ae3f8adf540889a) which made use of the line

`export("sol_stokes_3d_io", Th, sol, orderAndName, comm);`

and it looks like that fixed the issue for them. In the [most recent version of the example](https://github.com/FreeFem/FreeFem-sources/blob/master/examples/hpddm/stokes-io-3d.edp#L92), that line has been again revised to now read:

`savevtk("stokes-io-3d.vtu", Th, [u, uB, uC], uD, order = fforder, dataname = "velocity pressure");`

Is it advisable to switch back to the earlier `export` command to create a single VTU file as the output, or is there a way to alter my current version of the code to do so?

[stokes-3d-PETScagain.edp](https://community.freefem.org/uploads/short-url/Ao5JQT5ewCCvceiHrNkwk2kIAeP.edp) (2.0 KB)

---

<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:** [August 7, 2020, 2:41pm UTC](https://community.freefem.org/t/global-solution-with-petsc/239/9 "2020-08-07T14:41:41Z")

</div>

Do not use `export`. Use `savevtk` on the process which gathers the solution, e.g., `mpirank == 0`, with the additional parameter `communicator = mpiCommSelf`.

---

<div class="post-metadata">

**Author:** ![jgonzalez49](https://avatars.discourse-cdn.com/v4/letter/j/f05b48/32.png) [@jgonzalez49](https://community.freefem.org/u/jgonzalez49)\
**Post date:** [August 7, 2020, 3:21pm UTC](https://community.freefem.org/t/global-solution-with-petsc/239/10 "2020-08-07T15:21:48Z")

</div>

Thank you, passing that parameter to savevtk, while making sure it works on the `mpirank==0` process, does build a single mesh. However, it seems that while it contains the full geometry it only keeps the solution set from one of the four regions (at least from what I see in Paraview). Here is the current version of the script (I had to also alter the way the scalar and vector fields to match the [Stokes output example](https://github.com/FreeFem/FreeFem-sources/blob/master/examples/hpddm/stokes-io-3d.edp)):

[stokes-3d-PETScagain.edp](https://community.freefem.org/uploads/short-url/foWkUccH6zB706lMSM8OExsyRf0.edp) (2.1 KB)

and here are the PVD and VTU files it outputs:

[output.zip](https://community.freefem.org/uploads/short-url/1KtM4symrjZbfIm5mxooMaBRHHK.zip) (133.7 KB)

---

<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:** [August 7, 2020, 3:36pm UTC](https://community.freefem.org/t/global-solution-with-petsc/239/11 "2020-08-07T15:36:25Z")

</div>

Well, of course… you don’t save the reduced solution. Same for the `medit` command…

---

<div class="post-metadata">

**Author:** ![jgonzalez49](https://avatars.discourse-cdn.com/v4/letter/j/f05b48/32.png) [@jgonzalez49](https://community.freefem.org/u/jgonzalez49)\
**Post date:** [August 7, 2020, 3:45pm UTC](https://community.freefem.org/t/global-solution-with-petsc/239/12 "2020-08-07T15:45:44Z")

</div>

That did the trick, thank you very much for all your help! I attach the current version of the code for posterity.

[stokes-3d-PETScagain.edp](https://community.freefem.org/uploads/short-url/2z2WmVzeF9o0ONUtrO58vUeigPi.edp) (2.1 KB)