Sample Navier-Stokes Solver

General Discussion
1

Dear Friends,
I want to write a sample Navier Stokes solver using Rubin-Graves type linearization technique as u^(n+1)=u^n+u. The following is my code. I have some difficulties in the starting boundary values and also in the linearized non-linear terms. Would you please share your ideas and solutions.
Best regards,
Selcuk Han AYDIN, PhD.
===================== The Code ==================
// Parameters
int nn = 4;
real ut=1e-30; //If I set ut=1 results becomoes 1+e30
// Mesh
mesh Th = square(nn, nn);
// Fespace
// fespace Uh(Th, P1b);
//Fespace
fespace Uh(Th, P2);
Uh ux, uy, vx, vy, uxn, uyn;
fespace Ph(Th, P1);
Ph p, q, pn;
uxn=0; uyn=0; pn=0;

varf prob0(uxn,vn) = on(3,uxn=ut);
real[int] res=prob0(0,Uh);
uxn=res;
cout << uxn<<“:”<<uyn<< endl;

//Macro
macro Gradient(u) [dx(u), dy(u)] //
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //
macro UgradV(ux,uy,vx,vy) [ [ux,uy]'[dx(vx),dy(vx)] , [ux,uy]'[dx(vy),dy(vy)] ]//

problem LinNS([ux,uy,p],[vx,vy,q]) =
int2d(Th)(Gradient(ux)'*Gradient(vx)+Gradient(uy)'Gradient(vy)
+Gradient(uxn)‘Gradient(vx)+Gradient(uyn)‘Gradient(vy)
+ UgradV(uxn,uyn, ux, uy)'[vx,vy] //Error in these terms
+ UgradV(ux,uy,uxn,uyn)’[vx,vy]
+ UgradV(uxn,uyn,uxn,uyn)’[vx,vy]
- (p+pn) * Divergence(vx, vy)
+ Divergence(ux,uy)*q + Divergence(uxn,uyn)q
+1e-8(p+pn)*q )
+on(1,2,3,4, ux=0., uy=0.);

2

What are the reasons of the errors in the following code.

  1. Why there are -nan in the first case
  2. Why there are e+29 in the second case

Best regards;

mesh Th = square(4, 4);
fespace Vh(Th, P1);
Vh u, v;
func ue = (x^2+y^2-1)/4; //ue: exact solution

u=0;
solve p0(u,v) = on(1,2,3,4,u=ue)+on(0,u=0);

cout << “Solution u=” << u << endl;

varf prob0(u,v) = on(1,2,3,4,u=ue);
real[int] res=prob0(0,Vh);
u=res;
cout << “u=”<<u<< endl;