Segmentation Violation using PETSc

General Discussion
1

Dear FreeFem developers and users,
I’m writing you since I’m having some troubles once I am required to pass from local to global solution in my parallel script. I have an algorithm that is fine in sequential and is apparently working also in parallel, but once I want to move from the local to global solution the following error appears:

[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range

I’m copying here also my script that I’m developing trying to follow what explained in this lecture

load "msh3"
load "mat_edgeP1"
load "Element_Mixte3d"
include "macro_ddm.idp"
lockOrientation = false;
mesh3 Th("team4gmsh.mesh");
lockOrientation = true;

mesh3 ThPlt = Th;

int brickPlt = ThPlt(0.05,0.,0.).region;
plot(Th,cmm = "mesh from gmsh", wait = 1);
load "PETSc"
macro dimension()3// EOM            // 2D or 3D

buildDmesh(Th) //now the mesh is local

int brick = Th(0.05,0.,0.).region;
int outside = Th(0.7,0.0,0.0).region;

cout<<"Brick label = "<<brick<<endl;
cout<<"Outside label = "<<outside<<endl;

real mu = 4*pi*1.0e-7;
real nu=1.0/mu;
real rhoAl = 3.94e-8; //Ohm*m
real sigmaAl = 1.0/rhoAl;
func Pk = Edge03d;

fespace vhp(Th,P0);
fespace vh(Th,Pk);
vh [A0x,A0y,A0z] = [-0.05*y,0.05*x,0.];
vhp SSigma = 1+(sigmaAl-1)*(region == brick);

macro Curl(Ax,Ay,Az) [dy(Az)-dz(Ay), dz(Ax)-dx(Az), dx(Ay)-dy(Ax)]//EOM

real tempo;
real tau = 0.0119; //[s]
real dt = 1.0e-4; //[s];
real tend = 0.02; //[s];
int Imax=tend/dt;

func fact = exp(-tempo/tau);
int[int] lati = [17,18,19,20,21,22];

varf lhs([AX,AY,AZ],[Wx,Wy,Wz]) = int3d(Th)(nu*Curl(AX,AY,AZ)'*Curl(Wx,Wy,Wz))+int3d(Th)(SSigma*[Wx,Wy,Wz]'*[AX,AY,AZ]/dt)+on(lati,AX=-0.05*y*fact,AY=0.05*x*fact,AZ=0.);
varf rhs([AX,AY,AZ],[Wx,Wy,Wz]) = int3d(Th)(SSigma*[Wx,Wy,Wz]'*[A0x,A0y,A0z]/dt)+on(lati,AX=-0.05*y*fact,AY=0.05*x*fact,AZ=0.);

Mat AA;
int[int] n2o;
macro ThN2O()n2o// EOM
macro def(i)[i, i#B, i#C]// EOM
macro init(i)[i, i, i]// EOM
createMat(Th, AA, Pk)

fespace VhPlt(ThPlt, Pk); //LINE GENERATING PROBLEMS
int[int] SubIdx = restrict(vh,VhPlt,n2o); //LINE GENERATING PROBLEMS

real[int] rr(vh.ndof);

ofstream fpj("MPI_pj_TEAM4.dat");

for (int ii=1;ii<=Imax;ii++)
{
  real tempo = dt*ii;

  AA = lhs(vh,vh);
  set(AA, sparams = "-pc_type lu");
  rr = rhs(0,vh);
  vh [Ax,Ay,Az];
  Ax[] = AA^-1*rr;
VhPlt [Agx,Agy,Agz];
Agx[](SubIdx) = Ax[]; //LINE GENERATING PROBLEMS
[A0x,A0y,A0z] = [Ax,Ay,Az];
}

Here attached also the mesh I’m using.
https://drive.google.com/drive/folders/1zW7m6ir9tQzcvi8_EWUTAVQa8MI7X4l2?usp=sharing

Is someone ablo to understand where the problem is and why introducing the lines marked with the comment LINE GENERATING PROBLEM the error appears?

I’m already thanking so much anyone will try to help me!

yours sincerely,
Marco

2

Both int[int] n2o; and macro ThN2O() n2o // should be defined before buildDmesh(...) because the array n2o will be obtained when building distributed mesh. In your code, the n2o is an null array, so the error occurs when the restrict function is called.

3

Oh many thanks! That was actually a silly mistake! Thanks again!

By the way once corrected my script works perfectly fine since I’m setting -n 1, while the final outcome is totally wrong if I increase the number of processors and I actually do not know the reason. To be more precise increasing the number of processors the final results increases to (of some order of magnitudes!)

I have tried to follow exactly the path described in the tutorial by @prj, but something is still missing in my code such that the results is wrong.

Is someone able to find out the mistake?

(updated script here attached!)

load "msh3"
load "mat_edgeP1"
load "Element_Mixte3d"
include "macro_ddm.idp"

lockOrientation = false;
mesh3 Th("team4gmsh.mesh");
lockOrientation = true;

mesh3 ThPlt = Th;

int brickPlt = ThPlt(0.05,0.,0.).region;
plot(Th,cmm = "mesh from gmsh", wait = 1);

load "PETSc"
macro dimension()3// EOM            // 2D or 3D

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

buildDmesh(Th) //now the mesh is local

cout<<"Mesh labels = "<<labels(Th)<<endl;
cout<<"Mesh regions = "<<regions(Th)<<endl;

int brick = Th(0.05,0.,0.).region;
int outside = Th(0.7,0.0,0.0).region;

cout<<"Brick label = "<<brick<<endl;
cout<<"Outside label = "<<outside<<endl;

real mu = 4*pi*1.0e-7;
real nu=1.0/mu;
real rhoAl = 3.94e-8; //Ohm*m
real sigmaAl = 1.0/rhoAl;
func Pk = Edge03d;

fespace vhp(Th,P0);
fespace Vhpp(ThPlt,P0);
fespace vh(Th,Pk);
vh [A0x,A0y,A0z] = [-0.05*y,0.05*x,0.];
vhp SSigma = 1+(sigmaAl-1)*(region == brick);
Vhpp Sigma = 1+(sigmaAl-1)*(region == brickPlt);

macro Curl(Ax,Ay,Az) [dy(Az)-dz(Ay), dz(Ax)-dx(Az), dx(Ay)-dy(Ax)]//EOM

real tempo;
real tau = 0.0119; //[s]
real dt = 1.0e-4; //[s];
real tend = 0.02; //[s];
int Imax=tend/dt;

real fact;
int[int] lati = [17,18,19,20,21,22];

varf lhs([Ax,Ay,Az],[Wx,Wy,Wz]) = int3d(Th)(nu*Curl(Ax,Ay,Az)'*Curl(Wx,Wy,Wz))+int3d(Th)(SSigma*[Wx,Wy,Wz]'*[Ax,Ay,Az]/dt)+on(lati,Ax=-0.05*y*fact,Ay=0.05*x*fact,Az=0.);
varf rhs([Ax,Ay,Az],[Wx,Wy,Wz]) = int3d(Th)(SSigma*[Wx,Wy,Wz]'*[A0x,A0y,A0z]/dt)+on(lati,Ax=-0.05*y*fact,Ay=0.05*x*fact,Az=0.);

Mat AA;
macro def(i)[i, i#B, i#C]// EOM
macro init(i)[i, i, i]// EOM

createMat(Th, AA, Pk)



fespace VhPlt(ThPlt, Pk);

real[int] rr(vh.ndof);

ofstream fpj("MPI_pj_TEAM4.dat");

int i4 = 0.004/dt;
int i8 = 0.008/dt;
int i12 = 0.012/dt;
int i16 = 0.016/dt;
int i20 = 0.02/dt;
ofstream f4("MPIprof_4ms_gmsh.dat");
ofstream f8("MPIprof_8ms_gmsh.dat");
ofstream f12("MPIprof_12ms_gmsh.dat");
ofstream f16("MPIprof_16ms_gmsh.dat");
ofstream f20("MPIprof_20ms_gmsh.dat");
real step = 1e-3;;


for (int ii=1;ii<=Imax;ii++)
{
  real tempo = dt*ii;

  fact = exp(-tempo/tau);

  AA = lhs(vh,vh);
  set(AA, sparams = "-pc_type lu");
  rr = rhs(0,vh);
  vh [Ax,Ay,Az];
  Ax[] = AA^-1*rr;

real[int] AxPETSc;
  real[int] A0xPETSc;
  ChangeNumbering(AA,Ax[],AxPETSc);
  ChangeNumbering(AA,Ax[],AxPETSc,inverse = true);
  ChangeNumbering(AA,A0x[],A0xPETSc);
  ChangeNumbering(AA,A0x[],A0xPETSc,inverse = true);

  VhPlt [Agx,Agy,Agz];
  VhPlt [Ag0x,Ag0y,Ag0z];


  VhPlt [AgxS,AgyS,AgzS];
  VhPlt [Ag0xS,Ag0yS,Ag0zS];

Agx[](SubIdx) = Ax[];
  Ag0x[](SubIdx) = A0x[];

  mpiAllReduce(Agx[],AgxS[], mpiCommWorld, mpiSUM);
  mpiAllReduce(Ag0x[],Ag0xS[],mpiCommWorld,mpiSUM);

  VhPlt [Ex,Ey,Ez] = [(Ag0xS-AgxS)/dt,(Ag0yS-AgyS)/dt,(Ag0zS-AgzS)/dt];

  VhPlt [Jex,Jey,Jez] = Sigma*[Ex,Ey,Ez];

  real Pjoule = int3d(ThPlt,brickPlt)(1./Sigma*(Jex^2+Jey^2+Jez^2));

  VhPlt [Bx,By,Bz] = Curl(Agx,Agy,Agz);

  if (mpirank == 0)
  {
    {
      ofstream fpj("MPI_pj_TEAM4.dat",append);
      fpj<<tempo<<" "<<Pjoule<<endl;
    }



    if (ii==i4)
    {
      real zz = 0.0;
      while (zz<0.13)
      {
        {
          ofstream f4("MPIprof_4ms_gmsh.dat",append);
          f4<<zz<<" "<<Bz(0.,0.,zz)<<" "<<Agx(0.0,0.0,zz)<<" "<<Agy(0.0,0.0,zz)<<" "<<Agz(0.0,0.0,zz)<<endl;
        }
        zz=zz+step;
      }
    }

    if (ii==i8)
    {
      real zzz = 0.0;
      while (zzz<0.13)
      {
        {
          ofstream f8("MPIprof_8ms_gmsh.dat",append);
          f8<<zzz<<" "<<Bz(0.,0.,zzz)<<" "<<Agx(0.0,0.0,zzz)<<" "<<Agy(0.0,0.0,zzz)<<" "<<Agz(0.0,0.0,zzz)<<endl;
        }
        zzz=zzz+step;
      }
    }

    if (ii==i12)
    {
      real zzzz = 0.0;
      while (zzzz<0.13)
      {
        {
          ofstream f12("MPIprof_12ms_gmsh.dat",append);
          f12<<zzzz<<" "<<Bz(0.,0.,zzzz)<<" "<<Agx(0.0,0.0,zzzz)<<" "<<Agy(0.0,0.0,zzzz)<<" "<<Agz(0.0,0.0,zzzz)<<endl;
        }
        zzzz=zzzz+step;
      }
    }

    if (ii==i16)
    {
      real xx = 0.0;
      while (xx<0.13)
      {
        {
          ofstream f16("MPIprof_16ms_gmsh.dat",append);
          f16<<xx<<" "<<Bz(0.,0.,xx)<<" "<<Agx(0.0,0.0,xx)<<" "<<Agy(0.0,0.0,xx)<<" "<<Agz(0.0,0.0,xx)<<endl;
        }
        xx=xx+step;
      }
    }

    if (ii==i20)
    {
      real xxx = 0.0;
      while (xxx<0.13)
      {
        {
          ofstream f20("MPIprof_20ms_gmsh.dat",append);
          f20<<xxx<<" "<<Bz(0.,0.,xxx)<<" "<<Agx(0.0,0.0,xxx)<<" "<<Agy(0.0,0.0,xxx)<<" "<<Agz(0.0,0.0,xxx)<<endl;
        }
        xxx=xxx+step;
      }
    }
  }

  [A0x,A0y,A0z] = [Ax,Ay,Az];
  cout<<"Solved at t = "<<tempo<<" [s]"<<endl;



}
4

It seems you are using edge elements. Could you please try to add this line as well before buildDmesh, please?

5

Thanks for the extremely prompt answer!
By the way if I simply add the line you are mentioning I’m getting the following error:

the array size must be 1 not 3
 the array size must be 1 not 3

  current line = 428 mpirank 1 / 2
 Error line number 428, in file  macro: partitionPrivate in C:\Program Files (x86)\FreeFem++\\idp\macro_ddm.idp, before  token ;
Invalide array size  for  vectorial fespace function
  current line = 428 mpirank 0 / 2
Compile error : Invalide array size  for  vectorial fespace function
        line number :428, ;
error Compile error : Invalide array size  for  vectorial fespace function
        line number :428, ;
 code = 1 mpirank: 0

It seems like I’m messing up with array size…

6

Maybe try [Edge03ds0, Edge03ds0, Edge03ds0] then?

7

Unfortunately in such a way the initial error reappears:

[1]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range

8

What is SubIdx in your script? Why is the variable n2o never used?

9

Oh well… I’m sorry that was actually only an error of copy and paste… in my script SubIdx is defined just after the command createMat in the following way:

int[int] SubIdx = restrict(vh,VhPlt,n2o);

Sorry again for the typo!

10

OK, I see another problem, why are you using those macro def(u) and macro init(u)? Please just stick to what is on the distribution here.

11

I have tried to stick as close as possible to what is on the distribution, but the problem remains, everything is fine using 1 process while moving to -n 2 (or even more) there is a complete mess up of the results (not only the “post processed” ones, but already the varf solution array.).

I’m ampologising already for all my questions, but I’m actually not able to find out the mistake.

Here attached the final version of the script.

Thanks again for your help!

TestMPI_TEAM4.edp (5.2 KB)

12

Please provide a smaller mesh which reproduces the problem as well.

13

Here it is. The mesh is extremely small, so for sure not at convergence, but however should reproduce the problem in a proper way!

team4gmsh_lite.mesh (470.0 KB)

Thanks again

14

I see multiple mistakes:

  • you are defining your integers brick and outside after the call to buildDmesh, which is wrong because that means that processes may have different values;
  • you are not using the reduced array AgxS for the postprocessing but rather Agx, which is not correct;
  • you are using Ax at the end of your loop, but its values have changed in the ChangeNumbering call prior, so this messes up the subsequent iterations.
1 Like
15

Thank you very much Professor! Now my script works very well with any number of processors!

Your help was extremely precious!

Yours Sincerely,

Marco De Bastiani