Hi,
i try to write a code who calculate the solution of the problem:
\partial_t u + u \partial_x u - \nu \partial^2_x u =1+x, \ x \in [0,1], t \in [0,1]
with initial condition:
u(x,0)= -4 \nu \dfrac{x+3/2}{ (x+3/2)^2+\nu } - (x+2)
and boundary conditions:
u(0,t)= -4 \nu \dfrac{1+1/2 e^t}{ (1+1/2 e^t)^2 +\nu }-(1+e^t),
u(1,t)= -4 \nu \dfrac{ 2+1/2 e^t }{ (2+1/2 e^t)^2 +\nu} -(2+e^t).
The exact solution of this problem is:
u(x,t)=-4 \nu \dfrac{1+x+1/2 e^t} { (1+x+1/2 e^t)^2+\nu } -(1+x+e^t)
We can take \nu=1..
The goal is to calculate the approximate solution with FF+ and compare it to the exact solution compring the graphs of the two functions. I try to write the Following code but there is an error and i Don’t find how correct it.
Thank you in advance to the help.
load "iovtk"
verbosity=1;
real nu=1.;
int n = 100;
real a=0, b=1;
mesh Th=square(n,1,[a+x*(b-a),y/10]);
//plot(Th,wait=false);
fespace Vh(Th,P2,periodic=[[1,x],[3,x]]);
Vh u,v;
u = 0;
real T=0;
func f = 1+x;
func real uex(real t) //The exact solution
{
return -4*nu* (1.+x+(1./2.)*exp(t))/ ( (1.+x+(1./2.)*exp(t))^2+nu) - (1.+x+exp(t));
};
func up= -4*nu* (x+(3./2.))/((x+(3./2.))^2+nu)-(x+2.);//Initial condition
func real g(real t) //Boundary condition u(x=0,t)
{
return -4*nu*(1.+(1./2.)*exp(t))/((1+(1./2.)*exp(-t))^2+nu) -(1+exp(t));
};
func real r(real t) //Boundary condition u(x=1,t)
{
return -4*nu*(2.+(1./2.)*exp(t))/((2.+(1./2.)*exp(t))^2+nu) -(2+exp(t));
};
Vh fh=f;
Vh gh=g;
Vh rh=r;
Vh uph=up;
real dt = 0.1;
Vh g0,g1;
Vh nith = 0.0; // for nonlinear iterations
problem Pb(u,v) = int2d(Th)( u*v/dt + nu*dx(u)*dx(v) ) - int2d(Th)( up*v/dt + f*v)
+ int2d(Th) (g0*v) + int2d(Th) (g1*u*v) + int2d(Th) (nith*dx(u)*v)+ on(2,u=g)+ on(4,u=r);
int k=0;
int noit = 10; // maximum of iterations
real errnlpb = 0; // error of iteration of Pb
real normuh = 0;
ofstream fout("uh.txt");
fout << "# t uh" << endl;
ofstream fouti("U.txt");
fouti << "# t U" << endl;
for(int it=0; it<1.; ++it)
{
T += dt;
// Nonlinear iterations for Pb solution
nith = u; // start nonlinear iterations with last solution for Pb
g0=u*dx(u);
g1=dx(u);
cout << " Solving nonlinear problem.\n";
for(k=1;k<noit; ++k) {
up[]=u[];
Pb;
normuh = sqrt(int1d(Th,1)(square(u))); // L^2 nor of the solution of Pb
errnlpb = sqrt(int1d(Th,1)(square(u - nith))); // L^2 norm of the error in nonlinear iteration of Pb
cout << " " << k << ": error = " << errnlpb << endl;
if(errnlpb <= 1E-5*normuh) // stopping criteria. If relative error
break; // is less than 1E-5 stop nonlinear iterations of Pb
nith = u; // prepare new iteration
g0=u*dx(u);
g1=dx(u);
}
// Check convergence for Pb
if(k >= noit) {
cout << "Nonlinear solver for Pb did not converged.\n";
cout << "Error = " << errnlpb << endl;
cout << "|| u|| = " << normuh << endl;
exit(1); // exit the program with error code 1
}
Vh uh=uex;
plot(u,wait=1,dim=1,fill=true,cmm=it,value=1);
plot(uh,wait=1,dim=1,fill=true,cmm=it,value=1);
fout.flush;
fouti.flush;
}
SNESSolve shouldn’t know the difference between y and t…
