Defining PETSc matrix for Navier-Stokes

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 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");

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 `.

unformatted code 2*pi

Thanks.

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

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.

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

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 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 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:

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;

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.