Bonjour
j’ai le code suivant et je souhaite calculer
eh(t)= max_{x,y} |uh(x,y,t)-wh(x,y,t)|
où uh et wh sont deux solutions de deux problèmes différents calculer dans le code. Le voici joint et je vous remercie par avance pour toute aide.
load "UMFPACK64"
defaulttoUMFPACK64();
load "iovtk"
verbosity=1;
func real qdelta (real t, real s, real w) {
real t1 = t - rint(t);
real s1 = s - rint(s);
real delta=1.0/sqrt(w*pi);
if ( square(s1)+square(t1) <= square(delta))
return w;
else
return 0.0;
};
func real q (real t, real s, real w) {
real t1 = t - rint(t);
real s1 = s - rint(s);
real delta=1.0/sqrt(w*pi);
if ( square(s1)+square(t1) <= square(delta))
return 1;
else
return 0.0;
};
func real f (real x)
{
return 0.25*sin(2*pi*x);
};
// initial values for the spatial domain
int N=200;
real eps=0.5, w=3.0;
real t=0.00;
real dt=0.1;
int ni=5; // number of values of h=eps^2 w
real[int] h(ni); // value of h at each step
real[int] error(ni); // error in norm H1 t each step
real[int] alpha(ni); // alpha =1/log(2). log(error(h))/log(error(h/2))
ofstream fout("error.txt");
fout << "# step h error alpha" << endl;
for (int i=0;i<=ni-1;++i) {
mesh Th=square(N,N,[x,y]);
fespace Vh(Th,P1,periodic=[[2,y],[4,y],[1,x],[3,x]]);
fespace Wh(Th,P2,periodic=[[2,y],[4,y],[1,x],[3,x]]);
Vh uh=0.0;
Vh wh=0.0;
Vh vh;
Vh oldU=0.0;
Vh oldW=0.0;
real T=1.0;// arbitraire
int[int] Order= [1];
string DataName;
Vh fh=f;
Vh nith = 0.0; // for nonlinear iterations
Vh g0,g1,g2,g3;
Wh uplot;
// Nonlinear part is calculated in previous iteration (nith) and becomes a part of the right hand side
problem het(uh,vh,init=t) =
int2d(Th)((dx(uh)*dx(vh)+dy(uh)*dy(vh))*dt)
+int2d(Th)(uh*vh)
-int2d(Th)(oldU*vh)
+int2d(Th)(qdelta(x/eps,y/eps,w)*g0*vh*dt) +int2d(Th)(qdelta(x/eps,y/eps,w)*g1*uh*vh*dt)
-int2d(Th)(fh*vh*dt)
;
problem hom(wh,vh,init=t) =
int2d(Th)((dx(wh)*dx(vh)+dy(wh)*dy(vh))*dt)
+int2d(Th)(wh*vh)
-int2d(Th)(oldW*vh)
+int2d(Th)(g2*vh*dt) +int2d(Th)(g3*wh*vh*dt)
-int2d(Th)(fh*vh*dt)
;
uplot=q(x/eps,y/eps,w);
savevtk("q_"+i+".vtu",Th,uplot,dataname="qdelta", order=Order);
cout << "mean value of qdelta: " << int2d(Th)(uplot);
real max = -1e100; //initialisation de max
ofstream fouti("max.txt");//créer le fichiers où seront écrites les valeurs
for (t=0;t<=T;t+=dt) {
cout << "Time = " << t << endl;
fh=f(t);
oldU=uh;
oldW=wh;
int k=0;
int noit = 10; // maximum no of nonlinear iterations
real errnlhet = 0; // error of nonlinear iteration of het
real normuh = 0;
real errnlhom =0; // error of nonlinear iteration of hom
real normwh=0;
// Nonlinear iterations for heterogeneous solution
nith = uh; // start nonlinear iterations with last solution for het
g0=uh*uh/(1.0+uh)^2;
g1=1.0/(1.0+uh)^2;
cout << " Solving nonlinear heterogeneous problem.\n";
for(k=1;k<noit; ++k) {
het;
normuh = sqrt(int2d(Th)(square(uh))); // L^2 nor of the solution of het
errnlhet = sqrt(int2d(Th)(square(uh - nith))); // L^2 norm of the error in nonlinear iteration of het
cout << " " << k << ": error = " << errnlhet << endl;
if(errnlhet <= 1E-5*normuh) // stopping criteria. If relative error
break; // is less than 1E-5 stop nonlinear iterations of het
nith = uh; // prepare new iteration
g0=uh*uh/(1.0+uh)^2;
g1=1.0/(1.0+uh)^2;
}
// Check convergence for het
if(k >= noit) {
cout << "Nonlinear solver for het did not converged.\n";
cout << "Error = " << errnlhet << endl;
cout << "|| uh || = " << normuh << endl;
exit(1); // exit the program with error code 1
}
nith= wh; // start nonlinear iteration with last solution for hom
g2=wh*wh/(1.0+wh)^2;
g3=1.0/(1.0+wh)^2;
cout << " Solving nonlinear homogeneous problem.\n";
for(k=1;k<noit; ++k) {
hom;
normwh = sqrt(int2d(Th)(square(wh))); // L^2 nor of the solution of hom
errnlhom = sqrt(int2d(Th)(square(wh - nith))); // L^2 norm of the error in nonlinear iteration of hom
cout << " " << k << ": error = " << errnlhom << endl;
if (errnlhom <= 1E-5*normwh) // stopping criteria. If relative error
break; // is less than 1E-5 stop nonlinear iterations of het
nith = wh; // prepare new iteration
g2=wh*wh/(1.0+wh)^2;
g3=1.0/(1.0+wh)^2;
}
// Check convergence for het
if(k >= noit) {
cout << "Nonlinear solver for hom did not converged.\n";
cout << "Error = " << errnlhom << endl;
cout << "|| wh || = " << normwh << endl;
exit(1); // exit the program with error code 1
}
Vh eh=abs(uh-wh); // Définition de la fonction eh pour laquelle on veut calculer son max en (x,y)
real umaxt = eh[].max; // max au temps temp
umax = max( max, umaxt);// max jusqu'au temps t
fouti << umaxt << " "<< umax << endl;
}
savevtk("u_"+i+".vtu",Th,uh,dataname="uh", order=Order);
savevtk("w_"+i+".vtu",Th,wh,dataname="vh", order=Order);
real errH1 = sqrt(int2d(Th)(square(uh - wh)) + int2d(Th)(square(dx(uh)-dx(wh)) + square(dy(uh)-dy(wh))));
cout << " ==== errH1 ( eps = " << eps <<", omega = " << w << ", eps^2*omega = " << eps^2*w << ") = "<< errH1 <<endl;
h(i)=eps^2*w;
error(i)=errH1;
alpha(i)=i>0?1.0/log(2.0)*log(error(i))/log(error(i-1)):1.0;
fout << i <<" "<< h(i) <<" "<< error(i) <<" "<< alpha(i) << endl;
fout.flush;
eps=eps/2;
w=2*w;
}
