Hi,
i write the Following code to calculate the solution of two problems
du/dt-Delta u+ qdelta(x/eps,y/eps,z/eps) u/(1+u)= 0.25sin(2pi*x), (x,y,z) \in (0,1)^3, t \in [0,1)
avec la condition initiale u(x,0)=0 et les conditions au ord périodique.
qdelta(t,s,c)= 1 si (x,y,z) \in B_{delta}+ \Z^3 et 0 sinon.
le rayon de la boule B_{delta}= (3./4*(1./(pi*w)))^1/3.
et l’autre problème est
dw/dt-Delta w+w/(1+w)=0.25sin(2pi*x)
avec condition au bord periodique et condition initiale w(x,0)=0.
So i try with Following code:
Is my code correct? Please. Beacause when i visualize uh and wh there are static! is not normal. And if i try eps=0.5 and w=3, i see the ball in qdelta but the boundary of ball isn’t regular
Kind regards
load "UMFPACK64"
defaulttoUMFPACK64();
load "iovtk"
load "msh3"
verbosity=1;
func real qdelta (real t, real s, real c, real w) {
real t1 = t - rint(t);
real s1 = s - rint(s);
real c1 = c -rint(c);
real delta= ((3./4.)*(1./pi*w))^1/3;
if ( square(s1)+square(t1)+square(c1) <= square(delta))
return w;
else
return 0.0;
};
func real q (real t, real s, real c,real w) {
real t1 = t - rint(t);
real s1 = s - rint(s);
real c1 = c -rint(c);
real delta= ((3./4.)*(1./pi*w))^1/3;
if ( square(s1)+square(t1)+square(c1) <= 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=20;
real eps=0.1, w=2000;
real t=0.00;
real dt=0.1;
real[int] h; // value of h
real[int] error; // error in norm H1
ofstream fout("error.txt");
fout << "# h error" << endl;
ofstream fouti("errH1.txt");
fouti << "# t errH1" << endl;
mesh3 Th=cube(20,20,20);
fespace Vh(Th,P1,periodic=[[1,x,z],[3,x,z],[2,y,z],[4,y,z],[5,x,y],[6,x,y]]);
fespace Wh(Th,P2,periodic=[[1,x,z],[3,x,z],[2,y,z],[4,y,z],[5,x,y],[6,x,y]]);
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];
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) =
int3d(Th)((dx(uh)*dx(vh)+dy(uh)*dy(vh)+dz(uh)*dz(vh))*dt)
+int3d(Th)(uh*vh)
-int3d(Th)(oldU*vh)
+int3d(Th)(qdelta(x/eps,y/eps,z/eps,w)*g0*vh*dt) +int3d(Th)(qdelta(x/eps,y/eps,z/eps,w)*g1*uh*vh*dt)
-int3d(Th)(fh*vh*dt)
;
problem hom(wh,vh,init=t) =
int3d(Th)((dx(wh)*dx(vh)+dy(wh)*dy(vh)+dz(wh)*dz(vh))*dt)
+int3d(Th)(wh*vh)
-int3d(Th)(oldW*vh)
+int3d(Th)(g2*vh*dt) +int3d(Th)(g3*wh*vh*dt)
-int3d(Th)(fh*vh*dt)
;
uplot=q(x/eps,y/eps,z/eps,w);
//plot(uplot,fill=true,value=true,ps="q".eps");
savevtk("q.vtu",Th,uplot,dataname="qdelta", order=Order);
cout << "mean value of qdelta: " << int3d(Th)(uplot);
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(int3d(Th)(square(uh))); // L^2 nor of the solution of het
errnlhet = sqrt(int3d(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(int3d(Th)(square(wh))); // L^2 nor of the solution of hom
errnlhom = sqrt(int3d(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
}
real errH1 = sqrt(int3d(Th)(square(uh - wh)) + int3d(Th)(square(dx(uh)-dx(wh)) + square(dy(uh)-dy(wh))+ square(dz(uh)-dz(wh))));
cout << " ==== errH1 ( t = " << t <<", errH1 = " << errH1 << ") = "<< errH1 <<endl;
fout << t <<" "<< errH1 <<endl;
}
savevtk("u.vtu",Th,uh,dataname="uh", order=Order);
savevtk("w.vtu",Th,wh,dataname="vh", order=Order);
real errH1 = sqrt(int3d(Th)(square(uh - wh)) + int3d(Th)(square(dx(uh)-dx(wh)) + square(dy(uh)-dy(wh))+ square(dz(uh)-dz(wh))));
cout << " ==== errH1 ( eps = " << eps <<", omega = " << w << ", eps^2*omega = " << eps^2*w << ") = "<< errH1 <<endl;
h=eps^2*w;