///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef 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 of the License, or
/// (at your option) any later version.
///
/// Rheolef 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 Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================
//! @examplefile navier_stokes_solve.icc The Navier-Stokes equations with the method of characteristics -- solver function
using namespace std;
int navier_stokes_solve (
    Float Re, Float delta_t, field l0h, field& uh, field& ph,
    size_t& max_iter, Float& tol, odiststream *p_derr=0) {
  const space& Xh = uh.get_space();
  const space& Qh = ph.get_space();
  string label = "navier-stokes-" + Xh.get_geo().name();
  integrate_option iopt;
  iopt.set_family(integrate_option::gauss_lobatto);
  iopt.set_order(Xh.degree());
  trial u (Xh), p (Qh);
  test  v (Xh), q (Qh);
  form m  = integrate (dot(u,v), iopt);
  form a  = integrate (2*ddot(D(u),D(v)) + 1.5*(Re/delta_t)*dot(u,v), iopt);
  form b  = integrate (-div(u)*q, iopt);
  problem_mixed stokes (a, b);
  if (p_derr != 0) *p_derr << "[" << label << "] #n |du/dt|" << endl;
  field uh1 = uh;
  for (size_t n = 0; true; n++) { 
    field uh2 = uh1;
    uh1  = uh;
    field uh_star = 2.0*uh1 - uh2;
    characteristic X1 (    -delta_t*uh_star);
    characteristic X2 (-2.0*delta_t*uh_star);
    field l1h = integrate (dot(compose(uh1,X1),v), iopt);
    field l2h = integrate (dot(compose(uh2,X2),v), iopt);
    field lh  = l0h + (Re/delta_t)*(2*l1h - 0.5*l2h);
    stokes.solve (lh, field(Qh,0), uh, ph);
    field duh_dt = (3*uh - 4*uh1 + uh2)/(2*delta_t);
    Float residual = sqrt(m(duh_dt,duh_dt));
    if (p_derr != 0) *p_derr << "[" << label << "] "<< n << " " << residual << endl;
    if (residual < tol) {
      tol = residual;
      max_iter = n;
      return 0;
    }
    if (n == max_iter-1) {
      tol = residual;
      return 1;
    }
  }
}
