Dear All,
I pasted below and attached a 3d code for incompressible Navier-Stokes equations, discretized by a Crank-Nicolson type of discretization in time and the nonlinear term is treated by a fixed point iteration. The domain is the cude [-1,1]^3, and the problem that I am trying to simulate is the Eithier-Steinman flow problem. I do not have success in obtaining the first time step solution correctly, and after that the solver breaks down. I tried UMFPACK and GMRES. I would appreciate any help and share the successful code with anybode interested.
I am running successfully this type of problem, with the same discretization in 2d using freefem. I attached that code as well (Chorin_NSE.edp).
I thank you all in advance for your time and help,
Monika
Pasted code:
///////////////////////////////////////////////////////////////////////////////////////////
verbosity=5;
// loading
load "msh3"
load "tetgen"
load "medit"
//generation of domain=[-1,1]^3 and mesh/////////////////////////////////////////////////////////
real x0=-1.0,x1=1.0;
real y0=-1.0,y1=1.0;
int nn=5,mm=5;
mesh Thsq1=square(nn,mm,[x0+(x1-x0)*x,y0+(y1-y0)*y]);
func ZZ1min = -1.0;
func ZZ1max = 1.0;
func XX1 = x;
func YY1 = y;
int[int] ref31h = [0,1];
int[int] ref31b = [0,2];
mesh3 Th31h = movemesh23(Thsq1,transfo=[XX1,YY1,ZZ1max],refface=ref31h);
mesh3 Th31b = movemesh23(Thsq1,transfo=[XX1,YY1,ZZ1min],refface=ref31b);
func ZZ2 = x;
func YY2 = y;
func XX2min = -1.0;
func XX2max = 1.0;
int[int] ref32h = [0,3];
int[int] ref32b = [0,4];
mesh3 Th32h = movemesh23(Thsq1,transfo=[XX2min,YY2,ZZ2],refface=ref32h);
mesh3 Th32b = movemesh23(Thsq1,transfo=[XX2max,YY2,ZZ2],refface=ref32b);
func YY3min = -1.0;
func YY3max = 1.0;
func XX3 = x;
func ZZ3 = y;
int[int] ref33h = [0,5];
int[int] ref33b = [0,6];
mesh3 Th33h = movemesh23(Thsq1,transfo=[XX3,YY3max,ZZ3],refface=ref33h);
mesh3 Th33b = movemesh23(Thsq1,transfo=[XX3,YY3min,ZZ3],refface=ref33b);
////////////////////////////////
mesh3 Th = Th31h+Th31b+Th32h+Th32b+Th33h+Th33b; // "gluing" surface meshes to obtain the surface of cube
real[int] domaine = [0.5,0.5,0.5,145,0.001];
mesh3 Thfinal = tetg(Th,switch="pqaAAYYQ",nbofregions=1,regionlist=domaine); // Tetrahelize the interior of the cube with tetgen
//medit("tetg",Thfinal);
//exact solution for the Ethier-Steinman problem in 3d///////////////////////////////////
real a=1.25;
real d=1.0;
real nu=1.0;
func real u1e (real t) {
return -a*( exp(a*x)*sin(a*y+d*z) + exp(a*z)*cos(a*x+d*y) ) *exp(-nu*d*d*t);
}
func real u2e (real t) {
return -a*( exp(a*y)*sin(a*z+d*x) + exp(a*x)*cos(a*y+d*z) ) *exp(-nu*d*d*t);
}
func real u3e (real t) {
return -a*( exp(a*z)*sin(a*x+d*y) + exp(a*y)*cos(a*z+d*x) ) *exp(-nu*d*d*t);
}
func real pe (real t) {
return -0.5*(a*a)*( exp(2.0*a*x)*exp(2.0*a*y)*exp(2.0*a*z)+2.0*sin(a*x+d*y)*cos(a*z+d*x)*exp(a*(y+z))
+ 2.0*sin(a*y+d*z)*cos(a*x+d*y)*exp(a*(z+x))
+ 2.0*sin(a*z+d*x)*cos(a*y+d*z)*exp(a*(x+y)) )*exp(-2.0*nu*d*d*t);
}
/////////////////////////////////////////////////////////////////////////////////////
fespace Vh(Thfinal,P1b);
fespace Ph(Thfinal,P1);
Vh u1,u2,u3,v1,v2,v3;
Vh u1old,u2old,u3old,bd1,bd2,bd3;
Vh um1,um2,um3,e1,e2,e3;
Ph p,q;
u1=u1e(0); u2=u2e(0); u3=u3e(0); p=pe(0);
bd1=u1e(0); bd2=u2e(0); bd3=u3e(0);
u1old=0; u2old=0; u3old=0;
um1=0; um2=0; um3=0;
medit("solution", Thfinal, u1, u2,u3, wait=1);
real errLinfL2=0.;
real errL2=0., errGradL2=0.;
real err1, err2;
int itnl=0,iter=0, MaxNlIt=30;
real t=0;
real dt=0.005;
real alpha=1./dt;
int nbiter = 1; //200;
macro Grad(u1,u2,u3) [dx(u1),dx(u2),dx(u3),dy(u1),dy(u2),dy(u3),dz(u1),dz(u2),dz(u3)] //
macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) //
//Navier-Stokes with Dirichlet boundary conditions on on the sides of the box////////////////////////
macro NavierStokes(usesolver){
solve navierstokes([u1,u2,u3,p],[v1,v2,v3,q],solver=usesolver)=
int3d(Thfinal)( alpha*(u1*v1+u2*v2+u3*v3) )
+ int3d(Thfinal)( 0.5*nu*(Grad(u1,u2,u3)'*Grad(v1,v2,v3)) )
+ int3d(Thfinal)( 0.5*( u1old*dx(u1)*v1 + u2old*dy(u1)*v1 + u3old*dz(u1)*v1
+ u1old*dx(u2)*v2 + u2old*dy(u2)*v2 + u3old*dz(u2)*v2
+ u1old*dx(u3)*v3 + u2old*dy(u3)*v3 + u3old*dz(u3)*v3 ) )
- int3d(Thfinal)( div(v1,v2,v3)*p + div(u1,u2,u3)*q ) + int3d(Thfinal)( 1.e-5*p*q)
+int3d(Thfinal)( -alpha*(um1*v1+um2*v2+um3*v3) )
+ int3d(Thfinal)( 0.5*nu*(Grad(um1,um2,um3)'*Grad(v1,v2,v3)) )
+ int3d(Thfinal)( 0.5*(um1*dx(um1)*v1 + um2*dy(um1)*v1 + um3*dz(um1)*v1
+ um1*dx(um2)*v2 + um2*dy(um2)*v2 + um3*dz(um2)*v2
+ um1*dx(um3)*v3 + um2*dy(um3)*v3 + um3*dz(um3)*v3 ) )
+on(1,2,3,4,5,6,u1=bd1,u2=bd2,u3=bd3);
}//
//time loop///////////////////////////////////////////////////////////////////////
for (iter=1;iter<=nbiter;iter++)
{
t = t+dt;
um1 = u1;
um2 = u2;
um3 = u3;
bd1 = u1e(t);
bd2 = u2e(t);
bd3 = u3e(t);
err1 = 1;
cout << " iter = " << iter << " ------------------------ " << endl;
// fixed point iteration loop for nonlinear term ///////////////////////////////////
for (itnl=1; itnl<=MaxNlIt; itnl++)
{
cout<< "itnl = " << itnl << endl ;
u1old = u1;
u2old = u2;
u3old = u3;
NavierStokes(UMFPACK);
//NavierStokes(GMRES);
e1 = abs(u1 - u1old);
e2 = abs(u2 - u2old);
e3 = abs(u3 - u3old);
cout<< "e1[].max " << e1[].max << " e2[].max " << e2[].max<< " e3[].max " << e3[].max<<endl ;
cout<< "e1[].min " << e1[].min << " e2[].min " << e2[].min<< " e3[].min " << e3[].min<<endl ;
cout<< " " << endl ;
if ( (e1[].max < 0.0001) && (e2[].max < 0.0001) && (e3[].max < 0.0001))
{
break ;
}
medit("solution", Thfinal, u1, wait=1);
err1 = err1 + 1 ;
}
errL2 = int2d(Thfinal)((u1e(t)-u1)^2+(u2e(t)-u2)^2+(u3e(t)-u3)^2); //computes the L2 error
if ( errLinfL2 < errL2) {errLinfL2 = errL2;}
if ( err1 >= MaxNlIt ) { cout << "approx failed for iter = " << iter << endl; break; } //taking care of the fixed point iteration for the nonlinear loop
}
errLinfL2 = sqrt(errLinfL2);
cout << "errLinfL2= " << errLinfL2 << endl;
cout << "--------------------------------------" << endl;
medit("solution", Thfinal, u1, order=1);
///////////////////////////////////////////////////////////////////////////////////////////////////////////
--------------------------------------------------------------------
Monika Neda, Assistant Professor
Dept. of Mathematical Sciences
Univ. of Nevada Las Vegas (UNLV)
4505 Maryland PKWY, Box 454020
Las Vegas, NV 89154-4020
Email:
address@hiddenURL:
http://faculty.unlv.edu/neda/Phone: US +1 7028955170 Call
Fax: 702-895-4343
--------------------------------------------------------------------
Ethier_Steinman_NSE.edp