Singular matrix in PETSc

General Discussion
1

Hello,
can anyone help find out where is the problem here. See the code below :

matrix AA = [[A,B],
        [C,D]];
    //Mat DD = AA;
    //set(DD, sparams = "-pc_type lu -pc_factor_mat_solver_type mumps -mat_mumps_icntl_24 1 mat_mumps_icntl_33 -ksp_view");
    //the two commented lines are when the matrix is not invertible
    set(AA, sparams = "-pc_type lu -pc_factor_mat_solver_type mumps -mat_mumps_icntl_24 1 mat_mumps_icntl_33 -ksp_view");
    //here the matrix is invertible
    //notice that the only difference is the type
    Left=[Lefts, Leftpcap];
    xx=AA^-1 * Left;

I have the same matrix when I invert AA directly. The matrix is invertible and gives me the good solution, but when I define DD as Mat DD = AA and invert it, I get that the determinant is zero, and I don’t get any solution. I need the Mat type to use parallel computing in another code.

Thank you for your help,
Benfanich

2

No code to run => no help. Also, you are posting this in different posts: if you open another one, you will be timed out.

3

Sorry, about this, before posting this one I flagged the other one, if you are a mod please just remove it, because I didn’t describe my problem clearly, Here’s the code :

load "PETSc"
verbosity=0;
real cpu = clock();
int up = 1000, down = up -1, right = down-1, left = right-1;

////////////////////////Mesh: soil geometry////////////////////////////////////////////
//method parameters
real NN = 10, dt=0.1, t=0., tfin=100000, L=100; int i;
func real curve(real t) {
    return     L * (0.1 * (1 - cos(pi * (t/(L)))) + 0.45)/* (t <= 50) ? 50 : 50 + 0.5 * (t-50) */;
}

border DOWN(t=0,L){x=t; y=0; label=down;}
border right1(t=0, curve(L)){x=L; y=t; label = right;}
border curvee(t=L,0){x=t; y=curve(t);}	
border left1(t=curve(0),0){x=0; y=t; label=left;}
border left2(t=L,curve(0)){x=0; y=t; label = left;}
border UP(t=L,0){x=t; y=L; label=up;}
border right2(t=curve(L),L){x=L; y=t; label = right;}
real coef1=curve(L)*1./L, coef2=curve(0)*1./L, bla=1.0;
mesh Th = buildmesh(DOWN(NN)+right1(coef1*bla*NN)+curvee(bla*NN)+left1(coef2*bla*NN)+left2((1-coef2)*bla*NN)+UP(NN)+right2((1-coef1)*bla*NN));
fespace Ph(Th,P0);
Ph reg=region;
int bottom=reg(50,1), top=reg(50,99);

//soil parameters 
real q0, Ksup=0.25, Ksdown=2.00, Alphaup=0.028, Alphadown=0.016, hcapup=1/Alphaup, hcapdown= 1/Alphadown, thetasup=0.5, thetasdown=0.46, thetarup=0.12, thetardown=0.034
, nup=3, ndown=1.37
, mup=1-1/nup, mdown=1-1/ndown
, alphaup = 1/mup, alphadown= 1/mdown
, betaup= 1/nup, betadown= 1/ndown;
Ph condition1 = region == top , condition = 1-condition1/* region == bottom */ ; 
Ph Ks = Ksup*(condition1)   +   Ksdown*(condition); 
Ph Alpha = Alphaup*(condition1)   +   Alphadown*(condition);
Ph hcap = hcapup*(condition1)   +   hcapdown*(condition);
Ph thetar=thetarup*(condition1)   +   thetardown*(condition); 
Ph thetas=thetasup*(condition1)   +   thetasdown*(condition); 
Ph phi=thetas-thetar;
Ph n =nup*(condition1) + ndown*(condition);
Ph m = (1-1/nup)*(condition1) + (1-1/ndown)*(condition);
Ph alpha=alphaup*(condition1) + alphadown*(condition );
Ph beta=betaup*(condition1) + betadown*(condition);


macro KSS(s,t,r) (r*pow(abs(s),0.5)*pow(1-pow(abs(1-pow(abs(s),1./t)),t),2))//
macro KS(s) ((0<s && s<1) ? KSS(s,m,Ks) : (s>=1 ? Ks : 0)) //


macro norm(v) (sqrt(pow(dx(v),2)+1))//
macro grad(v) ([dx(v),dy(v)])//



///////////////////////////FE space////////////////////////////////////////////////
fespace Vh(Th, P1);
Vh s,v,h,w,sold, soldd, ls = (curve(x)-y), Delta = NN*(y==L/2.);
ls = ls/norm(ls);
fespace Wh(Th, [P0,P0]);
Wh [pcap1,pcap2],[zz1,zz2], [normal1,normal2];
[normal1,normal2] = grad(ls);
//////////////////////////Numerical Solution//////////////////////////////////////////

varf Amat(s,v)= int2d(Th, qft=qf1pTlump)(phi*s*v)+on(up,down,s=1);
varf Bmat([pcap1,pcap2],[v])=int2d(Th)(-(dt)*[pcap1,pcap2]'*grad(v));
varf Cmat([s], [zz1,zz2])=int2d(Th)((KS(sold))*hcap*grad(s)'*[zz1,zz2]);
varf Dmat([pcap1,pcap2],[zz1,zz2])=int2d(Th, qft=qf1pTlump)( (1e-12)* [pcap1,pcap2]'*[zz1,zz2]);

varf Leftsidepcap([pcap1,pcap2], [zz1,zz2]) = -int2d(Th)( (1e-12)* KS(sold)*zz2);
varf Leftsides(s,v) = int2d(Th, qft=qf1pTlump)(phi*sold*v) +on(up,down,s=1);



real[int] viso(41);
int ndofw=Wh.ndof, ndofv=Vh.ndof;
matrix A(ndofv, ndofv), B(ndofw,ndofv), C(ndofv, ndofw), D(ndofw,ndofw), AA(ndofw+ndofv,ndofw+ndofv);
real[int] Leftpcap(ndofw), Left(ndofv+ndofv),xx(ndofv + ndofw), Lefts(ndofv);
sold =1; 
for (int i = 0; i < viso.n; i++) viso[i] = 0+ (1-0)*(i*1./(viso.n-1));
Vh toplot=sold;
string vtkfolder=".//paraview_data/";
plot(toplot,wait=0,prev=1, dim=2, fill=1, value = 1,viso=viso(0:viso.n-1), nbiso=30,cmm=  "min: " +toplot[].min+" || max: "+toplot[].max+ "|| Time: "+t);
for(i=1; i*dt <= tfin ; i++){
    A=Amat(Vh,Vh);
    B=Bmat(Wh,Vh);
    C=Cmat(Vh,Wh);
    D=Dmat(Wh,Wh);
    Leftpcap = Leftsidepcap(0,Wh);
    Lefts = Leftsides(0,Vh);
    AA = [[A,B],
        [C,D]];
    //Mat DD = AA;
    //set(DD, sparams = "-pc_type lu -pc_factor_mat_solver_type mumps -mat_mumps_icntl_24 1 mat_mumps_icntl_33 -ksp_view");
    //the two commented lines are when the matrix is not invertible
    set(AA, sparams = "-pc_type lu -pc_factor_mat_solver_type mumps -mat_mumps_icntl_24 1 mat_mumps_icntl_33 -ksp_view");
    //here the matrix is invertible
    //notice that the only difference is the type
    Left=[Lefts, Leftpcap];
    xx=AA^-1 * Left;
    //xx=DD^-1 * Left;
    [s[], pcap1[]]=xx;
    t=i*dt;
    //richardPcap;
    sold=s;
    toplot=s;
    cout << "time : " << t << endl;
    if(i%1==0){plot(toplot,wait=0,prev=1, dim=2, fill=1, value = 1,viso=viso(0:viso.n-1), nbiso=30,cmm=  "min: " +toplot[].min+" || max: "+toplot[].max+ "|| Time: "+t);}
 } 
cpu = clock() - cpu;
cout << "Total time : " << cpu <<endl;

Check the comments in the for loop,

Thank you,

Benfanich

4

Here is the MUMPS error:

 On return from DMUMPS, INFOG(1)=              -9
 On return from DMUMPS, INFOG(2)=            7532

Try to run with -mat_mumps_icntl_14 1000.

5

It did not work but when I tried the value, 100000 it worked, but it’s very slow. Is there any way to speed it up? Thank you.

6

Use parallel computing or switch to a preconditioned iterative method instead of an exact factorization.