About the functions LinearCG

General Discussion
1

Hi everyone!
Just a small question about the precise nomenclature of the function LinearCG.
If I try to solve a very basic Poisson problem with homogeneous Dirichlet BC, the following piece of code works just fine:

varf a(u,v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
              + on(1,u=0);

varf rhs(u,v) = int2d(Th)(f()*v);

matrix A = a(Vh,Vh,solver=CG);
Vh b,u;
b[] = rhs(0,Vh);
u[] = A^-1*b[];

On the contrary, when I now use the routine LinearCG to achieve the same goal, it does not converge:

varf a(u,v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
              + on(1,u=0);

varf rhs(u,v) = int2d(Th)(f()*v);

matrix A = a(Vh,Vh);
Vh b,pp,u;
b[] = rhs(0,Vh);

/* Using conjugate gradient */
func real[int] opA(real[int] &p) {
  pp[] = A*p;
  return(pp[]);
}
u = 0;
LinearCG(opA,u[],b[],eps=1.e-6,nbiter=1000);

Yet, if I am not mistaken, the two operations should amount to the same, right?

Sorry in advance if there is something I misunderstood, and thank you for your help!

1 Like
2

Hi,

There is no misunderstanding. The only thing that you are missing is the contribution of boundary condition in your second example. Indeed, the reason of divergence of the solution is lack of boundary condition. You can fix this either by modifying your conjugate gradient function and turning it to a more advanced one (you may find a good sample on the Fluid-Structure Interaction page) or simply include the BC in your right hand side equation (varf rhs(u,v) = int2d(Th)(f()*v) + on(1,u=0); )

I personally prefer to combine rhs and lhs terms together to prevent such a problem (and for other reasons as well, e.g. keeping all BCs in one place). As a result, a working example would be like this, which produces the same results as if you use the simple solver=CG argument:

mesh Th = square(10,10);
fespace Vh(Th, P1);
func f = 5;

varf a(u,v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
              + int2d(Th)(f*v) + on(1,u=0);

matrix A = a(Vh,Vh);
Vh b,pp,u;
b[] = a(0,Vh);

func real[int] opA(real[int] &p) {
  pp[] = A*p;
  return(pp[]);
}

u = 0;
LinearCG(opA,u[],b[],eps=1.e-6,nbiter=100000);

plot(u, fill=true);
3

Dear Mojtaba,
thank you very much for your fast and very helpful answer.
I understand totally what you are saying… nevertheless, I am still a little puzzled.
Since the boundary condition I intend to impose is 0 boundary condition, nevertheless, I would think that whatever entry is on the right-hand side on this Dirichlet boundary does not matter too much…

In addition, running the example you propose, the result is “LinearCG does not converge”,
and the solution obtained, however not absurd, significantly differs from that obtained by entering simply:

u[] = A^-1*b[];

A more serious problem is reached if in the above example, I take

  • on(1,u=1)

as a Dirichlet boundary condition on the part 1 of the boundary. Then the LinearCG algorithm stops immediately.

Thank you again for your very valuable help in any event!!

4

Hi Charles,

I think part of the confusion is related to the application of iterative conjugate gradient methods and the situations for which these methods are used.

In addition, running the example you propose, the result is “LinearCG does not converge”,

Really??? Do you run the same code? For me, it produces a correct and equal result to the direct solver. Weird. Do you notice that I increased the number of iterations?

Regarding the boundary condition, you are totally right, and it would be a serious problem when you change your BC to a non-zero one (like “on(1,u=1)”). To deal with that, we may modify the code to be like this:

mesh Th = square(40,40);
fespace Vh(Th, P1);
func f = 5;

varf a(u,v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
              + int2d(Th)(f*v) + on(1,u=1);

matrix A = a(Vh,Vh,tgv=-1,solver=GMRES);
Vh u,pp,b;

b[] = a(0,Vh,tgv=-1);

func real[int] opA(real[int] &p) {
  pp[] = A*p;
  return pp[];
};

u[] = b[];
LinearCG(opA,u[],b[],eps=1.e-8,nbiter=1000);

plot(u, fill=true);

Important modifications are changing the solver (GMRES instead of CG) and initialization of u using the right hand side term, which again, considers the contribution of boundary conditions. This code works smoothly for me, so let me know if it doesn’t for you.

You may find this link useful in this matter. It contains some information about solving Poisson equation with LinearCG. Although their code is comprehensive and relatively complex, we can summarize their approach to be something like this:

int numnode=40;
mesh Th = square(numnode,numnode);

fespace Uh(Th,P1);
Uh u;

func f = 5;

varf a(u,v)=int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v)) + int2d(Th)(f*v) + on(1,u=1);
varf BC(u,v)=on(1,u=1);
matrix A=a(Uh,Uh,tgv=-1,solver=GMRES);
real[int] BC1=BC(0,Uh,tgv=-1);

func real[int] DJ(real[int] &u)
{
  real[int] Au=A*u;
  return Au;
};

real[int] b=a(0,Uh,tgv=-1);

u[] = BC1 ? b : 0;
u[] = b;

LinearCG(DJ,u[],b,eps=1.e-8,nbiter=10000);

plot(u,value=1,fill=1,cmm="Discrete solution");

The nice part of this code is the “u = BC1 ? b : 0;” statement. Indeed, it automatically takes into account whether you are using a 0 BC or not, and initializes the iteration accordingly.

5

Hallo both;

It was very clear that the BC was missing from the RHS. I have tried the first modified example by @mojtaba.barzegari and indeed does not work in ffpp 4.1 under win10. In the latest example, even if you keep matrix A = a(Uh, Uh, tgv=-1, solver=CG) all still works (I guess this just assembles half the matrix or?); hence, I assume that the important part of the code is the modification tgv=-1. Indeed, keeping default tgv results in a parasitic solution for me too. In the end, A is SPD and I do not see a good reason for LinearCG failing to converge. Could you please briefly explain why we need this tgv modification (or just provide a pointer where it is clarified)? Is that just because with the default tgv the residual has order tgv=1e30 components?

Regards,
/Fotis

6

Hi Fotios,

You are right, adding tgv=-1 is also very important, and I missed to mention it. But indeed, I cannot run my codes without initialization of u to rhs as well as using GMRES solver. If we use matrix A = a(Uh, Uh, tgv=-1, solver=CG) it still works and converges, but the solution is totally wrong (as you said).

Regarding the effect of tgv=-1, I had seen it in several FreeFem examples (in iterative examples and HPC examples when it is crucial to consider BC as a separate equation) and in the link I provided in my second answer. In FreeFem documentation, it is described as:

if tgv<0 then we put to 0 all term of the line i in the matrix, except diagonal term aii=1, and b[i]=“(Πhg)[i]”.

but to be honest, from mathematical perspective, I have no clue in my mind to relate this explanation to the fact that adding this argument causes the iterative gradient solver to converge.

Regarding what you said about the failure to run my provided codes (Charles had the same problem as well), I have to say it is super strange to me because all the codes I posted here converged and produced proper results in both Windows and Linux version of FreeFem v3.61. As a result, it seems that this could be a problem with v4.1, at least to say under windows.

7

Dear @mojtaba.barzegari thanks for your reply and the reference to the doc, it completely escaped my attention. My problem is that if tgv<0, then each row is replaced as described in the doc and the resulting matrix is not symmetric anymore; which is a bit counter-intuitive, but can be easily handled within the LinearCG implementation. I guess this is the source of my confusion here. Anyhow, I guess we all understand a bit more on LinearCG now.

Regards

8

Dear all,
thank you for all your answers, and for your all your help. Sorry to be late to get back to you. I have been contacted out of the forum by Pierre Jolivet (http://jolivet.perso.enseeiht.fr), who kindly explained to me the main stakes of the LinearCG function, and why my codes were not running. With his authorization, here are the main lines of his explanations.

In the piece of code

varf a(u,v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
+ int2d(Th)(f()v)
+ on(1,u=1);
matrix A = a(Vh,Vh,solver=CG);
Vh b,u;
b[] = a(0,Vh,tgv=tgv);
u[] = A^-1b[];

it is implicitly understood that the CG solver is used for the calculation of u[], together with a Jacobi (diagonal) preconditioner. The latter is necessary to “kill” the 1e30 values on the diagonal of A and the rows of b corresponding to Dirichlet boundary conditions in doing the matrix-vector products featured in CG.

Then I tried exactly the same thing, with a Jacobi preconditioner:

real tgv = 1.e30;

varf a(u,v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
+ int2d(Th)(f()*v)
+ on(1,u=1);

matrix A = a(Vh,Vh,tgv=tgv);
Vh b,u;
b[] = a(0,Vh,tgv=tgv);

/* Using conjugate gradient /
func real[int] opA(real[int] &p) {
pp[] = Ap;
return(pp[]);
}

/* Preconditioner */
func real[int] opQ(real[int] &p) {
for (int i=0; i<p.n; i++) {
p[i] = p[i] / A(i,i);
}
return§;
}
u[] = 0;
LinearCG(opA,u[],b[],precon=opQ,eps=1.e-6,nbiter=1000);

and it still does not work properly, because the CG algorithm stops after one iteration! The reason is that the initial residual, which is used in the CG algorithm to set the tolerance threshold for convergence is calculated, based upon the provided values of A, u and b, without preconditioning. Again, because of the tgv values, the threshold has an absurd value. To fix this, it is desirable to initialize u no longer at 0 but at 0, together with the correct values at the Dirichlet boundary (and 0 elsewhere), so that no tgv appears in the initial residual Au-b.

Hence, the token

u[] = a(0,Vh,tgv=-1);

before calling LinearCG makes things work fine.
The tgv=-1 means for the vector assembly that only boundary conditions are set without the tgv. When it comes to the matrix assembly, tgv=-1 or tgv=-2 imply that nonsymmetric or symmetric elimination are used, as described in the following tutorial by P. Jolivet:

http://jolivet.perso.enseeiht.fr/FreeFem-tutorial/#pf14

Eventually, Pierre advocates rather the use of LinearGMRES together with a nonsymmetric elimination of the Dirichlet degrees of freedom instead of the use of LinearCG together with symmetric elimination or use of the “tgv trick”.

Thank you very much (and to Pierre) for your help, in any event. This really made me understand a lot of things about the concrete implementation of routines in FreeFem!

1 Like