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;
lambdas(i)=lambda;
varf charge([u1, u2], [t1, t2])=int2d(Th)(
10^i*div(t1, t2)
);
real[int] F = charge(0,Vh);
toplot(0,i)=abs(F.sum);
// 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);
toplot(1,i)=abs(F0.sum);
// 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);
toplot(2,i)=abs(F1.sum);
// 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);
toplot(3,i)=abs(F2.sum);
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