Hello Everyone,
I am not getting the plot for [u1, u2].
kindly, suggest me to fix it.
code has been given as:
/*
coupled equation is
u_x + v_y = 0
((u,v).grad(u,v))* + grad(p) -1/Re * ( lap(u, v) ) = 0
(u, v). ( grad(T) ) - lap( T) = 0
B. C.
over square domain,
: all velocity components are zero
: temperature value is zero at left and right edge
and sin(pi*x) at the bottom edge
*/
int nn = 20;
real err;
real tol = 1.e-8;
real Re = 2000;
//verbosity =0;
mesh Th = square(nn, nn );
macro grad(u) [ dx(u), dy(u)] //
macro div(u1, u2) (dx(u1) + dy(u2) ) //
macro Grad(u1, u2) [ grad(u1), grad(u2) ] //
macro Ugradv( u1, u2, v1) ([u1, u2]'*grad(v1) ) //
macro UgradV(u1, u2, v1, v2) [Ugradv(u1, u2, v1), Ugradv(u1, u2, v2) ] //
// finite element spaces
fespace Vh( Th, [P2, P2, P1b, P2]);
Vh [u1, u2, p, T]; // variables
Vh [du1, du2, dp, dT]; // incremenet
Vh [v1, v2, q, TT]; // test function
problem heat([du1, du2, dp, dT], [v1, v2, q, TT]) =
int2d(Th) ( div(du1, du2)*q
+ [v1, v2]'*UgradV( du1, du2, u1, u2)
+ [v1, v2]'*UgradV( u1, u2, du1, du2)
+ 1./Re*(Grad(du1, du2): Grad(v1, v2))
- dp* div(v1, v2)
+ Ugradv(du1, du2, T)*TT
+ Ugradv(u1, u2, dT)*TT
+ grad(dT)'*grad(TT)
- q* div(du1, du2)
- 1.e-8*dp*q
)
- int2d(Th) (
[v1, v2]'*UgradV(u1, u2, u1, u2)
+[v1, v2]'*UgradV(u1, u2, u1, u2)
+ 1./Re*( Grad(u1, u2): Grad(v1, v2) )
- p* div(v1, v1)
+ TT* Ugradv(u1, u2, T)
+ grad(T)'*grad(TT)
- q*div(u1, u2)
- 1.e-8*p*q
)
+ on( 1, 2,3,4, du1=0, du2 = 0)
+ on( 2, 4, dT = 0)
+ on(1, dT = sin(3.14*x));
[u1, u2, p, T] = [0,0,0,0];
for (int step = 0; step < 15; step++)
{
heat;
u1[] = u1[] - du1[];
real solerr = u1[].linfty;
real stperr = du1[].linfty;
err = stperr/solerr;
cout<< "step is == "<< step<<endl;
cout<< "error is == "<< err<<endl;
if(err< tol) break;
if( step > 3 && err > 1.e3) break;
plot([u1, u2],p,T, wait = 0, value =1 );
//plot( T, wait =0, fill =1, value =1);
}