# Sample Navier-Stokes Solver

**URL:** https://community.freefem.org/t/sample-navier-stokes-solver/2164
**Category:** General Discussion
**Created:** [December 4, 2022, 9:47am UTC](https://community.freefem.org/t/sample-navier-stokes-solver/2164 "2022-12-04T09:47:08Z")
**Posts on this page:** 2
**Page:** 1

<div class="post-metadata">

### Author: ![shaydin](https://avatars.discourse-cdn.com/v4/letter/s/bc79bd/32.png) [@shaydin](https://community.freefem.org/u/shaydin)
#### Post date: [December 4, 2022, 9:47am UTC](https://community.freefem.org/t/sample-navier-stokes-solver/2164/1 "2022-12-04T09:47:08Z")

</div>

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.);

---

<div class="post-metadata">

### Author: ![shaydin](https://avatars.discourse-cdn.com/v4/letter/s/bc79bd/32.png) [@shaydin](https://community.freefem.org/u/shaydin)
#### Post date: [December 4, 2022, 6:32pm UTC](https://community.freefem.org/t/sample-navier-stokes-solver/2164/2 "2022-12-04T18:32:54Z")

</div>

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;