// -*- c++ -*-

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2003  Laboratoire J.-L. Lions UPMC Paris

//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2, or (at your option)
//  any later version.

//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

//  $Id: navier-stokes.ff,v 1.5 2007/02/14 21:34:23 delpinux Exp $

// Navier-Stokes flow around a cylinder
// Uses a projection scheme over the velocity divergence free space.

// Fictitious domain method uses the poor elimination method to run
// faster in this example.

// Only 5 iterations are performed

vector n = (50,10,10);
vertex A = (0,0,0);
vertex B = (5,1,1);
mesh M = structured(n,A,B);
scene S = pov("navier-stokes.pov");
domain Omega = domain(S,outside(<1,0,0>));

mesh obstacle = surface(Omega,M);

femfunction Ki(M) = 1 - one(<1,0,0>);

// Computes the area of Omega
double area = int[M](Ki);

// Cylinder diameter is 0.2 => Re = 400
function uin = 12.5;
function mu = 0.005;
femfunction u0(M) = uin;
femfunction v0(M) = 0;
femfunction w0(M) = 0;
femfunction p0(M) = 2 - x;

double dt = 0.01;
double eps = 0.001;
double pp;
double qq;

double i = 1;

function u=0;
function v=0;
function w=0;

do{

    cout << "========== Navier-Stokes step " << i << "\n";

    solve(u) in Omega by M memory(matrix=none),method(type=eliminate),
    krylov(type=cg,precond=diagonal){

	pde(u)
	    (1/dt)*u - div(mu*grad(u))
	    = ((1/dt)*convect([u0,v0,w0],dt,u0) - dx(p0));
	u = 0 on obstacle;
	u = uin on M xmin;
	u = convect([u0,v0,w0],dt,u0) on M xmax;
    };

    solve(v) in Omega by M memory(matrix=none),method(type=eliminate),
    krylov(type=cg,precond= diagonal){
	
	pde(v)
	    (1/dt)*v - div(mu*grad(v))
	    = ((1/dt)*convect([u0,v0,w0],dt,v0) - dy(p0));
	v = 0 on obstacle;
	v = 0 on M xmin;
	v = convect([u0,v0,w0],dt,v0) on M xmax;
	v = 0 on M ymin;
	v = 0 on M ymax;
    };

    solve(w) in Omega by M memory(matrix=none),method(type=eliminate),
    krylov(type=cg,precond=diagonal){

	pde(w)
	    (1/dt)*w - div(mu*grad(w))
	    = ((1/dt)*convect([u0,v0,w0],dt,w0) - dz(p0));
	w = 0 on obstacle;
	w = 0 on M xmin;
	w = convect([u0,v0,w0],dt,w0) on M xmax;
	w = 0 on M zmin;
	w = 0 on M zmax;
    };

    qq = int[M](dx(u)+dy(v)+dz(w))/area;
    solve(q) in Omega by M memory(matrix=none),method(type=eliminate),
    krylov(type=cg,precond=diagonal){

	pde(q)
	    (eps*dt)*q - div(dt*grad(q)) = (dx(u)+dy(v)+dz(w) - qq);

	q = 0 on M xmax;
    };

    p0 = p0 - q;
    pp = int[M](p0)/area;
    p0 = p0 - pp;
    u0 = u + (dt*dx(q));
    v0 = v + (dt*dy(q));
    w0 = w + (dt*dz(q));

    i = i+1;

} while(i <= 5);

// Final results.

save(medit,"velocity",M);
save(medit,"velocity",[u,v,w],M);
save(medit,"pressure",M);
save(medit,"pressure",p0,M);