navier_stokes_solve.icc

Go to the documentation of this file.

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;
    }
  }
}