search

FreeFem / FreeFem-sources

FreeFEM source code
https://freefem.org/
Other
776 stars 191 forks source link

Memory leak in parallel MUMPS #175

Closed abaikin closed 3 years ago

abaikin commented 3 years ago

Dear FF developers,

I found that when I solve the SLAE using MUMPS with MPI in the time loop, it leads to a memory leak. My FF ver. is 4.8. In v. 3.62 all was OK.

load "MUMPS"
load "metis"
/* Domain geometry */
real ROut = 100.0; /* [m], Radius of outer circle */
real RIn = 0.1; /* [m], Radius of inner circle */

/* Boundary conditions for inner and outer circles */
real pOut = 23e6;
real pIn = 12e6;

/* Mesh definition */
int NN = 50*2*5*2;
int NOut = NN;
int NIn = NN/10;
int labelIn = 1; /* Label of outer circle */
int labelOut = 2; /* Label of inner circle */

/* Reservoir parameters */

real Ct = 1e-10; /* [1/Pa], Reservoir compressibility */
real kr = 1e-14; /* [m^2], Reservoir permeability */
real eta = 5e-3; /* [Pa*s], Oil viscosity */

real phi0  = 0.2; /* [], Porosity */

real dt = 20; /* [s], */
real tMax = 20000; /* [s], */

/* Inner circle border definition */
border GammaIn(t = 0.0,2.0*pi){ 

    x = RIn*cos(t);
    y = RIn*sin(t);
    label = labelIn;
};

/* Outer circle border definition */
border GammaOut(t = 0.0,2.0*pi){
    x = ROut*cos(t);
    y = ROut*sin(t);
    label = labelOut;
}; 

mesh Th = buildmesh( GammaOut(NOut) + GammaIn(-NIn),nbvx=100000 ); 

fespace Ph0(Th,P0);

Ph0 part;
if(mpisize > 1){

  metisdual(part[], Th, mpisize);
  broadcast(processor(0, mpiCommWorld), part[]);

}

Th = change(Th, fregion = part[][nuTriangle] );

/* Finite element space definition */
fespace Vh(Th, P1); // P2 for better convergence
Vh p, phi;
Vh pOld;
pOld = pOut;
/* Bilinear form definition for matrix */
varf FiltrationBilinearForm(p,phi) = 
    int2d(Th, mpirank)(
        phi0*Ct*p*phi
    )
    + int2d(Th,mpirank)(
        dt*kr/eta*(dx(p) * dx(phi) + dy(p) * dy(phi)) 
    ) 
    + on(labelIn, p = pIn) 
    + on(labelOut, p = pOut);

/* Linear form definition for rhs */
varf RhsLinearForm(p,phi) = 
    int2d(Th, mpirank)(
        phi0*Ct*pOld*phi
    )
    + on(labelIn, p = pIn)
    + on(labelOut, p = pOut);

real t = 0;
matrix A;
real [int] b(Vh.ndof); 
while( t < tMax ){

    A = FiltrationBilinearForm(Vh,Vh); /* Matrix creation */

    b = RhsLinearForm(0,Vh); /* Rhs creation */ 

    {
        set(A, solver = sparsesolver, master = -1);

        p[] = A^-1*b ; /* SLAE solution */
    }

    cout << "t = " << t << endl;

    pOld = p;
    t = t + dt;
}
prj- commented 3 years ago

My advice is to use MUMPS through the PETSc interface, which is not leaking, and more efficient.

load "PETSc"
include "macro_ddm.idp"
/* Domain geometry */
real ROut = 100.0; /* [m], Radius of outer circle */
real RIn = 0.1; /* [m], Radius of inner circle */

/* Boundary conditions for inner and outer circles */
real pOut = 23e6;
real pIn = 12e6;

/* Mesh definition */
int NN = 50*2*5*2;
int NOut = NN;
int NIn = NN/10;
int labelIn = 1; /* Label of outer circle */
int labelOut = 2; /* Label of inner circle */

/* Reservoir parameters */

real Ct = 1e-10; /* [1/Pa], Reservoir compressibility */
real kr = 1e-14; /* [m^2], Reservoir permeability */
real eta = 5e-3; /* [Pa*s], Oil viscosity */

real phi0  = 0.2; /* [], Porosity */

real dt = 20; /* [s], */
real tMax = 20000; /* [s], */

/* Inner circle border definition */
border GammaIn(t = 0.0,2.0*pi){

    x = RIn*cos(t);
    y = RIn*sin(t);
    label = labelIn;
};

/* Outer circle border definition */
border GammaOut(t = 0.0,2.0*pi){
    x = ROut*cos(t);
    y = ROut*sin(t);
    label = labelOut;
};

mesh Th = buildmesh( GammaOut(NOut) + GammaIn(-NIn),nbvx=100000 );

Mat A;
createMat(Th, A, P1);

/* Finite element space definition */
fespace Vh(Th, P1); // P2 for better convergence
Vh p, phi;
Vh pOld;
pOld = pOut;
/* Bilinear form definition for matrix */
varf FiltrationBilinearForm(p,phi) =
    int2d(Th)(
        phi0*Ct*p*phi
    )
    + int2d(Th)(
        dt*kr/eta*(dx(p) * dx(phi) + dy(p) * dy(phi))
    )
    + on(labelIn, p = pIn)
    + on(labelOut, p = pOut);

/* Linear form definition for rhs */
varf RhsLinearForm(p,phi) =
    int2d(Th)(
        phi0*Ct*pOld*phi
    )
    + on(labelIn, p = pIn)
    + on(labelOut, p = pOut);

real t = 0;
real [int] b(Vh.ndof);
set(A, sparams = "-pc_type lu");
while( t < tMax ){

    A = FiltrationBilinearForm(Vh,Vh); /* Matrix creation */

    b = RhsLinearForm(0,Vh); /* Rhs creation */

        p[] = A^-1*b ; /* SLAE solution */

    cout << "t = " << t << endl;
    plotD(Th, p, cmm = "Solution");
    pOld = p;
    t = t + dt;
}
frederichecht commented 3 years ago

Dear user,

effectively , we missing to clean 2 pointeur in this case ??? I will correct rapidly !

Best Regards,

Frédéric Hecht.

Laboratoire Jacques-Louis Lions, UPMC Sorbonne Université BC187, 4 Place Jussieu, 75252 PARIS cedex 05, France Campus Jussieu, Barre 15-25, 3 etage Bureau 307 Projet Alpines , Inria de Paris, 2 rue Simone Iff Voie DQ12 75012 Paris tel: +33 1 44274411, mob: +33 6 62198986, fax: +33 1 44277200 @.*** https://www.ljll.math.upmc.fr/hecht software: FreeFem++ web site: http://www.freefem.org/ff++

abaikin commented 3 years ago

Dear Pierre, thank you for the advice. I will prefer to use it for large time-consuming problems.

Dear Dr. Hecht, I will appreciate it if you fix this bug. Sometimes, it is more convenient to use old MUMPS interface than PetSc distributed solution.

frederichecht commented 3 years ago

I have put the correction in the develop branch

diff --git a/plugin/mpi/MUMPS.cpp b/plugin/mpi/MUMPS.cpp index b044b003..201056fb 100644 --- a/plugin/mpi/MUMPS.cpp +++ b/plugin/mpi/MUMPS.cpp @@ -122,6 +122,15 @@ public: delete [] id.irn; delete [] id.jcn; delete [] id.a; +

Le 28 avr. 2021 à 08:34, abaikin @.***> a écrit :

Dear FF developers,

I found that when I solve the SLAE using MUMPS with MPI it leads to a memory leak. My FF ver. is 4.8. In v. 3.62 all was OK.

load "MUMPS" load "metis" / Domain geometry / real ROut = 100.0; / [m], Radius of outer circle / real RIn = 0.1; / [m], Radius of inner circle /

/ Boundary conditions for inner and outer circles / real pOut = 23e6; real pIn = 12e6;

/ Mesh definition / int NN = 50252; int NOut = NN; int NIn = NN/10; int labelIn = 1; / Label of outer circle / int labelOut = 2; / Label of inner circle */

/ Reservoir parameters /

real Ct = 1e-10; / [1/Pa], Reservoir compressibility / real kr = 1e-14; / [m^2], Reservoir permeability / real eta = 5e-3; / [Pas], Oil viscosity */

real phi0 = 0.2; / [], Porosity /

real dt = 20; / [s], / real tMax = 20000; / [s], /

/ Inner circle border definition / border GammaIn(t = 0.0,2.0*pi){

x = RIncos(t); y = RInsin(t); label = labelIn; };

/ Outer circle border definition / border GammaOut(t = 0.0,2.0pi){ x = ROutcos(t); y = ROut*sin(t); label = labelOut; };

mesh Th = buildmesh( GammaOut(NOut) + GammaIn(-NIn),nbvx=100000 );

fespace Ph0(Th,P0);

Ph0 part; if(mpisize > 1){

metisdual(part[], Th, mpisize); broadcast(processor(0, mpiCommWorld), part[]);

}

Th = change(Th, fregion = part[][nuTriangle] );

/ Finite element space definition / fespace Vh(Th, P1); // P2 for better convergence Vh p, phi; Vh pOld; pOld = pOut; / Bilinear form definition for matrix / varf FiltrationBilinearForm(p,phi) = int2d(Th, mpirank)( phi0Ctp*phi )

  • int2d(Th,mpirank)( dtkr/eta(dx(p) dx(phi) + dy(p) dy(phi)) )
  • on(labelIn, p = pIn)
  • on(labelOut, p = pOut);

/ Linear form definition for rhs / varf RhsLinearForm(p,phi) = int2d(Th, mpirank)( phi0CtpOld*phi )

  • on(labelIn, p = pIn)
  • on(labelOut, p = pOut);

real t = 0; matrix A; real [int] b(Vh.ndof); while( t < tMax ){

A = FiltrationBilinearForm(Vh,Vh); / Matrix creation /

b = RhsLinearForm(0,Vh); / Rhs creation /

{ set(A, solver = sparsesolver, master = -1);

  p[] = A^-1*b ; /* SLAE solution */

}

cout << "t = " << t << endl;

pOld = p; t = t + dt; }

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/FreeFem/FreeFem-sources/issues/175, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABLPNDA4DWZWGKZCKNJVS33TK6T6HANCNFSM43WNACLQ.

abaikin commented 3 years ago

Thank you! Now it is OK.