# Defining PETSc matrix for Navier-Stokes

**URL:** <https://community.freefem.org/t/defining-petsc-matrix-for-navier-stokes/210>\
**Category:** General Discussion\
**Created:** [January 7, 2020, 3:40pm UTC](https://community.freefem.org/t/defining-petsc-matrix-for-navier-stokes/210 "2020-01-07T15:40:59Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![david9333](https://avatars.discourse-cdn.com/v4/letter/d/65b543/32.png) [@david9333](https://community.freefem.org/u/david9333)\
**Post date:** [January 7, 2020, 3:40pm UTC](https://community.freefem.org/t/defining-petsc-matrix-for-navier-stokes/210/1 "2020-01-07T15:40:59Z")

</div>

Hi!

I’m new to FreeFem and parallel computing. I attended the tutorial sessions of the FreeFem Days in December, when PETSc was extensively used.  
I am currently on a coupled Vlasov-Navier-Stokes system and I want numerical simulations. I began with solving Navier-Stokes. I have a running code for sequential computing but it is too slow. Therefore, I am trying to switch to parallel computing.

My code generates an error at the line where I am building the PETSc matrix. To sum up, I have Pk=[P2,P2,P2,P1] and a 3d mesh called cylindre. The line “buildMat(cylindre,getARGV(”-split",1),A,Pk,mpiCommWorld)" generates : “Error line number 241, in file macro: partition in /usr/local/ff++/mpich3/lib/ff++/4.4-3/idp/macro\_ddm.idp”

I found this syntax in the [stokes-fieldsplit-3d-PETSc.edp](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/stokes-fieldsplit-3d-PETSc.edp) example (which runs smoothly on my computer) and I do not understand what I did differently.

Thank you all for your help!

David

Here is the full (minimal) code.

```
macro dimension()3//
include "macro_ddm.idp"
load "msh3"
load "tetgen"
load "PETSc"

macro def(i)[i, i#B, i#C,i#D]// // vector field definition

//Mesh
real r00 = 9.00e-1;
real L00 = 9.40;
int m=2;
int nlayer = L00*m;

border c(t=0,2*pi){x=r00*cos(t);y=r00*sin(t);};
mesh cercle = buildmesh(c(2*pi*r00*m));

int[int] rup=[0,1], rdown=[0,2], rmid=[1,3];
mesh3 cylindre = buildlayers(cercle,nlayer,
  coef = 1.,
  zbound = [0,L00],
  labelmid = rmid,
  reffaceup = rup,
  reffacelow = rdown);

buildDmesh(cylindre)

// Physical data
real u0 = 0.5e2;
real rhoair = 1.18e-3;
real nu = 1e-5*1e4;
real eps = 1e-8;

real dt = 0.01;

//Finite element spaces
func Pk=[P2,P2,P2,P1];
fespace Ph(cylindre,P0);
fespace Uh(cylindre,P1);
fespace Vh(cylindre,[P2,P2,P2,P1]);

//Navier-Stokes
func u1entree = 0;
func u2entree = 0;
func u3entree = u0*(1-(x^2+y^2)/(r00^2)); // Poiseuille at entry of the cylinder

Vh [u1,u2,u3,p];
Vh [v1,v2,v3,ptest];

// Initial data
Vh [up1,up2,up3,pini];
up1=0;
up2=0;
up3 = u0*(1-(x^2+y^2)/(r00^2));

varf bilinNS([u1,u2,u3,p],[v1,v2,v3,ptest]) = int3d(cylindre)(
  rhoair*(u1*v1 + u2*v2 + u3*v3)/dt
      +nu*(dx(u1)*dx(v1)+dy(u1)*dy(v1)+dz(u1)*dz(v1)
            +dx(u2)*dx(v2)+dy(u2)*dy(v2)+dz(u2)*dz(v2)
            +dx(u3)*dx(v3)+dy(u3)*dy(v3)+dz(u3)*dz(v3)
            )
      - p*(dx(v1)+dy(v2)+dz(v3))
      - ptest*(dx(u1)+dy(u2)+dz(u3))
      + p*ptest*eps
  )+ 
int3d(cylindre)(rhoair*(
      -convect([up1,up2,up3],-dt,up1)*v1/dt
        -convect([up1,up2,up3],-dt,up2)*v2/dt
        -convect([up1,up2,up3],-dt,up3)*v3/dt
  ))+
on(1, u1=u1entree, u2=u2entree, u3=u3entree)+on(3, u1=0, u2=0, u3=0);

matrix<real> Loc = bilinNS(Vh,Vh);
real[int] rhs = bilinNS(0,Vh);

Mat A;
**buildMat(cylindre,getARGV("-split",1), A, Pk, mpiCommWorld)**
A = Loc;

//From here on I am copying the example stokes-fieldsplit-3d-PETSc.edp, I do not really know what is happening
Vh<real> def(sol) = [1.0,1.0,1.0,2.0];
string[int] names(2);
names[0] = "velocity";
names[1] = "pressure";

set(A, sparams = "-ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -pc_fieldsplit_detect_saddle_point -fieldsplit_velocity_sub_pc_type lu " + " -fieldsplit_pressure_sub_pc_type lu -fieldsplit_velocity_sub_pc_factor_mat_solver_type mumps " + " -fieldsplit_pressure_sub_pc_factor_mat_solver_type mumps -ksp_monitor -ksp_view " + " -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_ksp_max_it 5 -fieldsplit_pressure_ksp_type gmres -fieldsplit_pressure_ksp_max_it 5 -ksp_rtol 1e-6", fields = sol[], names = names)

sol[] = 0.0;
sol[] = A^-1 * rhs;

plotD(cylindre,sol,cmm="sol distribuee");
```

---

<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:** [January 7, 2020, 5:36pm UTC](https://community.freefem.org/t/defining-petsc-matrix-for-navier-stokes/210/2 "2020-01-07T17:36:43Z")

</div>

Could you please paste your code in a non-formatted way? It’s impossible to copy/paste it right now. Surround the code by three `.

```auto
unformatted code 2*pi

```

Thanks.

---

<div class="post-metadata">

**Author:** ![david9333](https://avatars.discourse-cdn.com/v4/letter/d/65b543/32.png) [@david9333](https://community.freefem.org/u/david9333)\
**Post date:** [January 7, 2020, 9:48pm UTC](https://community.freefem.org/t/defining-petsc-matrix-for-navier-stokes/210/3 "2020-01-07T21:48:38Z")

</div>

Sorry about that, I’ve edited my previous message!  
Thanks for your answer!

---

<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:** [January 7, 2020, 9:55pm UTC](https://community.freefem.org/t/defining-petsc-matrix-for-navier-stokes/210/4 "2020-01-07T21:55:53Z")

</div>

The macro `init` is missing.  
Right after `macro def`, add something like  
`macro init(u)[u, u, u, u]//`  
There is a semicolon missing at the end of the `set(A, sparams = "...")` line as well.  
The other errors seem to come from your side, for example, it’s not correct to write the following.

```auto
Vh [up1,up2,up3,pini];
up1=0;
up2=0;
up3 = u0*(1-(x^2+y^2)/(r00^2));

```

---

<div class="post-metadata">

**Author:** ![david9333](https://avatars.discourse-cdn.com/v4/letter/d/65b543/32.png) [@david9333](https://community.freefem.org/u/david9333)\
**Post date:** [January 8, 2020, 4:23pm UTC](https://community.freefem.org/t/defining-petsc-matrix-for-navier-stokes/210/5 "2020-01-08T16:23:50Z")

</div>

Thank you for your answer!  
I’ve added the macro init and fixed the mistakes you pointed out. The command `ff-mpirun -np 1 file_name.edp -fglut name -wg` does not produce an error anymore.  
But when I change 1 into 2, I have the following errors (repeated 3 times), and then I have to kill the process (nothing happens anymore, but the execution isn’t finished). Could you please explain what it means?  
Thank you again so much!

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------  
[0]PETSC ERROR: Arguments are incompatible  
[0]PETSC ERROR: Local row size 1734 cannot be larger than global row size 768  
[0]PETSC ERROR: See [https://www.mcs.anl.gov/petsc/documentation/faq.html](https://www.mcs.anl.gov/petsc/documentation/faq.html) for trouble shooting.  
[0]PETSC ERROR: Petsc Release Version 3.12.2, Nov, 22, 2019  
[0]PETSC ERROR: /usr/local/ff++/mpich3/bin/FreeFem+±mpi on a named [pmichel.ljll.math.upmc.fr](http://pmichel.ljll.math.upmc.fr) by davidmichel Wed Jan 8 16:54:11 2020  
[0]PETSC ERROR: Configure options MAKEFLAGS= --prefix=/usr/local/ff++/mpich3/ff-petsc//real --with-debugging=0 --with-cc=/usr/local/mpich3/bin/mpicc --with-cxx=/usr/local/mpich3/bin/mpic++ --with-fc=/usr/local/mpich3/bin/mpif90 --with-cxx-dialect=C++11 COPTFLAGS="-O3 -mtune=generic" CXXOPTFLAGS="-O3 -mtune=generic" FOPTFLAGS="-O3 -mtune=generic" --with-ssl=0 --with-x=0 --with-fortran-bindings=0 --with-scalar-type=real --with-blaslapack-include= --with-blaslapack-lib="-Wl,-rpath,/opt/intel/mkl/lib -L/opt/intel/mkl/lib -lmkl\_rt -lmkl\_sequential -lmkl\_core -liomp5 -lpthread -lmkl\_rt -lmkl\_sequential -lmkl\_core -liomp5 -lpthread -lmkl\_rt -lmkl\_sequential -lmkl\_core -liomp5 -lpthread -lmkl\_rt -lmkl\_sequential -lmkl\_core -liomp5 -lpthread" --download-scalapack --download-metis --download-ptscotch --download-mumps --download-hypre --download-parmetis --download-superlu --download-suitesparse --download-tetgen --download-hpddm --download-hpddm-commit=83c462a --download-slepc PETSC\_ARCH=fr  
[0]PETSC ERROR: #1 MatSetSizes() line 154 in /Users/hecht/ff/ff-git/3rdparty/ff-petsc/petsc-3.12.2/src/mat/utils/gcreate.c

Here is the code I used:

```auto
include "macro_ddm.idp"
load "msh3"
load "tetgen"
load "PETSc"

macro def(i)[i, i#B, i#C,i#D]// // Définition d'un champ de vecteur
macro init(u)[u,u,u,u]//

//Mesh
real r00 = 9.00e-1;
real L00 = 9.40;
int m=2;
int nlayer = L00*m;

border c(t=0,2*pi){x=r00*cos(t);y=r00*sin(t);};
mesh cercle = buildmesh(c(2*pi*r00*m));

int[int] rup=[0,1], rdown=[0,2], rmid=[1,3];
mesh3 cylindre = buildlayers(cercle,nlayer,
  coef = 1.,
  zbound = [0,L00],
  labelmid = rmid,
  reffaceup = rup,
  reffacelow = rdown);

buildDmesh(cylindre)

// Physical data
real u0 = 0.5e2;
real rhoair = 1.18e-3;
real nu = 1e-5*1e4;
real eps = 1e-8;

real dt = 0.01;

//Finite element spaces
func Pk=[P2,P2,P2,P1];
fespace Ph(cylindre,P0);
fespace Uh(cylindre,P1);
fespace Vh(cylindre,[P2,P2,P2,P1]);

func u1entree = 0;
func u2entree = 0;
func u3entree = u0*(1-(x^2+y^2)/(r00^2)); // Poiseuille at entry of the cylinder

Vh [u1,u2,u3,p];
Vh [v1,v2,v3,ptest];

// Initial data
Vh [up1,up2,up3,pini] = [0,0,u0*(1-(x^2+y^2)/(r00^2)),0];

varf bilinNS([u1,u2,u3,p],[v1,v2,v3,ptest]) = int3d(cylindre)(
  rhoair*(u1*v1 + u2*v2 + u3*v3)/dt
      +nu*(dx(u1)*dx(v1)+dy(u1)*dy(v1)+dz(u1)*dz(v1)
            +dx(u2)*dx(v2)+dy(u2)*dy(v2)+dz(u2)*dz(v2)
            +dx(u3)*dx(v3)+dy(u3)*dy(v3)+dz(u3)*dz(v3)
            )
      - p*(dx(v1)+dy(v2)+dz(v3))
      - ptest*(dx(u1)+dy(u2)+dz(u3))
      + p*ptest*eps
  )
  + int3d(cylindre)(rhoair*(
      -convect([up1,up2,up3],-dt,up1)*v1/dt
        -convect([up1,up2,up3],-dt,up2)*v2/dt
        -convect([up1,up2,up3],-dt,up3)*v3/dt
  ))
  + on(1, u1=u1entree, u2=u2entree, u3=u3entree)
  + on(3, u1=0, u2=0, u3=0);

matrix<real> Loc = bilinNS(Vh,Vh);

real[int] rhs = bilinNS(0,Vh);

Mat A;
buildMat(cylindre,getARGV("-split",1), A, Pk, mpiCommWorld)

A = Loc;

//From here on I am copying the example stokes-fieldsplit-3d-PETSc.edp, I do not really know what is happening
Vh<real> def(sol) = [1.0,1.0,1.0,2.0];
string[int] names(2);
names[0] = "vitesse";
names[1] = "pression";

set(A, sparams = "-ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -pc_fieldsplit_detect_saddle_point -fieldsplit_velocity_sub_pc_type lu " + " -fieldsplit_pressure_sub_pc_type lu -fieldsplit_velocity_sub_pc_factor_mat_solver_type mumps " + " -fieldsplit_pressure_sub_pc_factor_mat_solver_type mumps -ksp_monitor -ksp_view " + " -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_ksp_max_it 5 -fieldsplit_pressure_ksp_type gmres -fieldsplit_pressure_ksp_max_it 5 -ksp_rtol 1e-6", fields = sol[], names = names);

sol[] = 0.0;
sol[] = A^-1 * rhs;

```

---

<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:** [January 8, 2020, 5:22pm UTC](https://community.freefem.org/t/defining-petsc-matrix-for-navier-stokes/210/6 "2020-01-08T17:22:08Z")

</div>

So, a couple of small wrong things:

1. you are missing a `macro dimension()3//` before `include "macro_ddm.idp"`, by default the value is 2.
2. `buildMat` is a combination of `buildDmesh` and `createMat`. So, you either want to replace `buildMat(cylindre,getARGV("-split",1), A, Pk, mpiCommWorld)` by `createMat(cylindre, A, Pk)`, or remove `buildDmesh` and call `buildMat(cylindre,getARGV("-split",1), A, Pk, mpiCommWorld)` only (before you do the finite element assemblies, e.g., after the definition of `func Pk`)
3. if you don’t know what the `-pc_type fieldsplit` do, why use it? Why don’t use a plain `-pc_type lu`?

Let me know if you need some more help.

---

<div class="post-metadata">

**Author:** ![Qianrui](https://avatars.discourse-cdn.com/v4/letter/q/9fc348/32.png) [@Qianrui](https://community.freefem.org/u/Qianrui)\
**Post date:** [May 23, 2022, 2:02am UTC](https://community.freefem.org/t/defining-petsc-matrix-for-navier-stokes/210/7 "2022-05-23T02:02:47Z")

</div>

> [@prj](#):
>
> getARGV(“-split”,1),

> getARGV(“-split”,1),

What does this parameter do here?