Compute Residuals

General Discussion
1

Hi everyone,

I’m trying to estimate the residual error when solving a pde with solve. More precisely freefem is solving a system which looks like : find U st AU = B where A is a matrix, B a vector. I would like to compute |AU-B|/|B|. I come up with the following code, but the results are far from the one i expected. Could someone explain to me what I’ve been doing wrong ? I suspect that i misunderstood how to use varf…

Thanks a lot,

Simon


real h = 0.01;
int Nm=floor(1./h); //Number of FE nodes on a unit segment
mesh Th=square(Nm,Nm,[x,y]);
fespace Vh(Th,P1);
Vh u,v,f;
f=1;

{
solve VA(u,v, solver=sparsesolver)=
       int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v))
   - int2d(Th)(f*v) 
   + on(1, 2, 3, 4, u=0);
}

varf bilinear(u,v) = int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v))
    	+ on(1, 2, 3, 4, u=0);

varf linear(unused,v) = int2d(Th)(f*v);

matrix A = bilinear(Vh, Vh);
Vh U, B;
U[] = A*u[];
B[] = linear(0, Vh);
U[] -= B[];	
cout << "residuals (expected small): " << ", " << sqrt(u[]'*u[])/sqrt(B[]'*B[]) << endl;
2

Its not a good idea to mix case on vars that look the same lol,

real h = 0.01;
int Nm=floor(1./h); //Number of FE nodes on a unit segment
mesh Th=square(Nm,Nm,[x,y]);
fespace Vh(Th,P1);
Vh u,v,f;
f=1;

{
solve VA(u,v, solver=sparsesolver)=
       int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v))
   - int2d(Th)(f*v) 
   + on(1, 2, 3, 4, u=0);
}

varf bilinear(u,v) = int2d(Th)(dx(u)*dx(v) + dy(u)*dy(v))
    	+ on(1, 2, 3, 4, u=0);

varf linear(unused,v) = int2d(Th)(f*v);

matrix A = bilinear(Vh, Vh);
Vh U, B;
U[] = A*u[];
B[] = linear(0, Vh);
U[] -= B[];	
cout << "residuals (expected small): " << ", " << sqrt(u[]'*u[])/sqrt(B[]'*B[]) << endl;
cout << "residuals (expected small): " << ", " << sqrt(U[]'*U[])/sqrt(B[]'*B[]) << endl;
plot(u,cmm="lc",wait=1,fill=1);
plot(U,cmm="ic",wait=1,fill=1);



type or paste code here
3

Indeed :slight_smile: thanks for the help !
The issue remains though : for this code example, the relative error is about 10% (which i find quite big, but i may be wrong).

4

You are missing + on(1, 2, 3, 4, unused=0) in your linear variational formulation.

5

Okk, yes now it works ! I did know that i was supposed to specify it on the linear part.
Thanks a lot :slight_smile:

6

Well, I guess in general a mistake I made early on is that a dot product is not
an integral. You would porlbably be better off doing something like `

real top=int2d(Th)(UU);
real bottom=int2d(Th)(BB);
cout<<" top “<<top<<” bottom “<<bottom<<” l2 "<<sqrt(top/bottom)<<endl;
`

7

I see but i’m looking for the error stemming from the resolution of the matricial system AU=B, where U is the unknown and represents the coefficients in the decomposition of the solution of the PDE. I see U as a big vector of R^n which contains the residual terms a(u,\phi_i)-b(\phi_i) for all i (vertex of the mesh) and want to sum (in a L2 sense) these local errors to get a “macroscopic” error.
This error is obviously linked to the one you are mentioning, but i guess this is not the one I’m looking for.