← Back to branch summary

~ubuntu-branches/ubuntu/trusty/rheolef/trusty

« back to all changes in this revision

Viewing changes to doc/pexamples/navier_stokes_solve.icc

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito
  • Date: 2012-04-06 09:12:21 UTC
  • mfrom: (1.1.5)
  • Revision ID: package-import@ubuntu.com-20120406091221-m58me99p1nxqui49
Tags: 6.0-1
http://bugs.debian.org/667356
http://bugs.debian.org/661689
* New upstream release 6.0 (major changes):
  - massively distributed and parallel support
  - full FEM characteristic method (Lagrange-Gakerkin method) support
  - enhanced users documentation 
  - source code supports g++-4.7 (closes: #667356)
* debian/control: dependencies for MPI distributed solvers added
* debian/rules: build commands simplified
* debian/librheolef-dev.install: man1/* to man9/* added
* debian/changelog: package description rewritted (closes: #661689)

Show diffs side-by-side

added added

removed removed

Lines of Context:
doc/pexamples/navier_stokes_solve.icc
 
1
using namespace std;
 
2
int navier_stokes_solve (
 
3
    Float Re, Float delta_t, field l0h, field& uh, field& ph,
 
4
    size_t& max_iter, Float& tol, odiststream *p_derr=0) {
 
5
  const space& Xh = uh.get_space();
 
6
  const space& Qh = ph.get_space();
 
7
  string label = "navier-stokes-" + Xh.get_geo().name();
 
8
  quadrature_option_type qopt;
 
9
  qopt.set_family(quadrature_option_type::gauss_lobatto);
 
10
  qopt.set_order(Xh.degree());
 
11
  form m (Xh, Xh, "mass", qopt);
 
12
  form a (Xh, Xh, "2D_D");
 
13
  form mp(Qh, Qh, "mass");
 
14
  a = a + 1.5*(Re/delta_t)*m;
 
15
  solver sa (a.uu());
 
16
  form b (Xh, Qh, "div"); b = -b;
 
17
  solver_abtb stokes (a.uu(), b.uu(), mp.uu());
 
18
  if (p_derr != 0) *p_derr << "[" << label << "] #n |du/dt|" << endl;
 
19
  field uh1 = uh;
 
20
  for (size_t n = 0; true; n++) { 
 
21
    field uh2 = uh1;
 
22
    uh1  = uh;
 
23
    field uh_star = 2.0*uh1 - uh2;
 
24
    characteristic X1 (    -delta_t*uh_star);
 
25
    characteristic X2 (-2.0*delta_t*uh_star);
 
26
    field l1h = riesz(Xh, compose(uh1,X1), qopt);
 
27
    field l2h = riesz(Xh, compose(uh2,X2), qopt);
 
28
    field lh  = l0h + (Re/delta_t)*(2*l1h - 0.5*l2h);
 
29
    stokes.solve (lh.u() - a.ub()*uh.b(), -(b.ub()*uh.b()),
 
30
                  uh.set_u(), ph.set_u());
 
31
    field duh_dt = (3*uh - 4*uh1 + uh2)/(2*delta_t);
 
32
    Float residual = sqrt(m(duh_dt,duh_dt));
 
33
    if (p_derr != 0) *p_derr << "[" << label << "] "<< n << " " << residual << endl;
 
34
    if (residual < tol) {
 
35
      tol = residual;
 
36
      max_iter = n;
 
37
      return 0;
 
38
    }
 
39
    if (n == max_iter-1) {
 
40
      tol = residual;
 
41
      return 1;
 
42
    }
 
43
  }
 
44
}