Unexpected values by int1d depending on datatype of a scalar multiplier

Dear community

I am really at a loss concerning the behavior of int1d when applying a scalar multiplier to the integral.

I was comparing two models in 2d linear elasticity which where mathematically supposed to yield the same results, but didn’t.
I tracked the error down to the integration of the linear part of the problem, which contains a material parameter lambda.
In one example it is defined as

fespace Wh(Th, P2,periodic=[[2, y], [4, y], [1, x], [3, x]]);
Wh lambda=10;

and in the other as

real lambda = 10;

Contrary to what I expected, this small change of datatype lead to a great difference in the result of my linear part of variational problem, defined as:

varf Macro([u1,u2],[v1,v2]) = int2d(Th)(lambda*(div(v1,v2)));
real[int] F = MacroL111(0,Vh);

So i did a small study, where I set varying powers of 10 for lambda and 4 different datatypes (real,P0 Fe-function, P1 Fe-function and P2 Fe function). The program is the following

// Mesh
mesh Th = square(10, 10);

// Fespace for solution
fespace Vh(Th, [P2,P2],periodic=[[2, y], [4, y], [1, x], [3, x]]);
Vh [t1, t2];
// Fespaces to hold material parameters
fespace Wh0(Th, P0,periodic=[[2, y], [4, y], [1, x], [3, x]]);
fespace Wh1(Th, P1,periodic=[[2, y], [4, y], [1, x], [3, x]]);
fespace Wh2(Th, P2,periodic=[[2, y], [4, y], [1, x], [3, x]]);
// how many powers of lambda to test
int lim=20;
real[int,int] toplot(4,lim);
real[int] lambdas(lim);
// Macro
macro div(u,v) ( dx(u)+dy(v) ) //
ofstream file("data.txt");
for (int i=0;i<lim;i++){

	// using a constant value
	real lambda=10^i;
	varf charge([u1, u2], [t1, t2])=int2d(Th)(
			 10^i*div(t1, t2)
	real[int] F = charge(0,Vh);

	// using a P0-function
	Wh0 lambda0=lambda;
	varf charge0([u1, u2], [t1, t2])=int2d(Th)(
			 lambda0*div(t1, t2)
	real[int] F0 = charge0(0,Vh);

	// using a P1-function
	Wh1 lambda1=lambda;
	varf charge1([u1, u2], [t1, t2])=int2d(Th)(
			 lambda1*div(t1, t2)
	real[int] F1 = charge1(0,Vh);

	// using a P2-function
	Wh2 lambda2=lambda;
	varf charge2([u1, u2], [t1, t2])=int2d(Th)(
			 lambda2*div(t1, t2)
	real[int] F2 = charge2(0,Vh);
	cout << "i="<<i<<"\n";
	cout << "lambda,lambda0,lambda1,lambda2="<<toplot(0,i)<<" "<<toplot(1,i)<<" "<<toplot(2,i)<<" "<<toplot(3,i)<<"\n";
	file << lambda<<" "<< toplot(0,i)<<" "<<toplot(1,i)<<" "<<toplot(2,i)<<" "<<toplot(3,i)<<"\n";

Normally, i would expect each row of the written “data.txt”(attached below) to have the same values (allowing for a bit of numerical noise) and each column to have values increasing by a factor of 10. Neither of those is the case however.

What is the right datatype to give to material parameters? real or Fe-function?
And why am i not getting a factor 10 between the entries in a given column in “data.txt” regardless of datatype used.

Thanks in advance for your help!

Here are my values data.txt:
1 6.11766e-16 6.11766e-16 6.11766e-16 7.34931e-16
10 4.97144e-15 4.97144e-15 4.97144e-15 1.37977e-14
100 3.91589e-14 3.91589e-14 3.91589e-14 4.20455e-14
1000 1.13519e-12 1.13519e-12 6.05835e-13 1.16761e-12
10000 2.82316e-12 2.82316e-12 2.82316e-12 6.06324e-12
100000 2.58975e-11 2.58975e-11 2.58975e-11 5.8412e-11
1e+06 5.75486e-10 5.75486e-10 1.64394e-10 1.05953e-10
1e+07 5.75168e-09 5.75168e-09 5.75168e-09 1.00041e-09
1e+08 2.81437e-08 2.81437e-08 2.81437e-08 3.80681e-08
1e+09 9.49917e-08 9.49917e-08 3.93015e-07 4.8432e-08
1e+10 7.07837e-07 7.07837e-07 7.07837e-07 1.27408e-06
1e+11 9.22311e-05 9.22311e-05 9.57368e-06 9.22311e-05
1e+12 0.000399005 0.000399005 0.000517261 0.000927341
1e+13 0.00716078 0.00716078 0.00240004 0.00199759
1e+14 0.00967624 0.00967624 0.00711276 0.00961087
1e+15 0.228921 0.228921 0.142984 0.211343
1e+16 3.33969 3.33969 2.10562 1.21469
1e+17 25.3452 25.3452 25.3452 78.3452
1e+18 520.327 520.327 838.327 1218.33
-9.22337e+18 5642.55 5642.55 5642.55 6778.55

It seems that you are computing
where (t1,t2) is perdiodic. Theoretically you shoud get 0.
What you observe is just random roundoff error (that you mutiply by a zooming factor).

1 Like

It is not linked to the fact of it being periodic (i tested).
But you are right in that the integral should be zero (one can apply the divergence theorem) and i am just scaling numerical noise!
If instead of the sum of my vector i output the max i am getting reasonable values in all cases!
Thanks for the tip!

You can see this in 1D- the integral of the derivative is just the difference in
values at the end which is probably zero in your case.
If its periodic you can expand it as fourier and see the DC term drops out
leaving a zero integral.
I would think numerical
noise would be easy to spot but that may not always be the case if
the problem or solution has issues.

Are you on Window?

Try to use 10.^i instead of 10^i.