rheolef  6.5
navier_stokes_cavity.cc
Go to the documentation of this file.
1 #include "rheolef.h"
2 using namespace rheolef;
3 using namespace std;
4 #include "navier_stokes_solve.icc"
5 #include "navier_stokes_criterion.icc"
6 #include "cavity.icc"
7 int main (int argc, char**argv) {
8  environment rheolef (argc, argv);
9  if (argc < 2) {
10  cerr << "usage: " << argv[0] << " <geo> <Re> <err> <n_adapt>" << endl;
11  exit (1);
12  }
13  geo omega (argv[1]);
14  adapt_option_type options;
15  Float Re = (argc > 2) ? atof(argv[2]) : 100;
16  options.err = (argc > 3) ? atof(argv[3]) : 1e-2;
17  size_t n_adapt = (argc > 4) ? atoi(argv[4]) : 5;
18  Float delta_t = 0.05;
19  options.hmin = 0.004;
20  options.hmax = 0.1;
21  space Xh = cavity_space (omega, "P2");
22  space Qh (omega, "P1");
23  field uh = cavity_field (Xh, 1.0);
24  field ph (Qh, 0);
25  field fh (Xh, 0);
26  for (size_t i = 0; true; i++) {
27  size_t max_iter = 1000;
28  Float tol = 1e-5;
29  navier_stokes_solve (Re, delta_t, fh, uh, ph, max_iter, tol, &derr);
30  odiststream o (omega.name(), "field");
31  o << catchmark("Re") << Re << endl
32  << catchmark("delta_t") << delta_t << endl
33  << catchmark("u") << uh
34  << catchmark("p") << ph;
35  o.close();
36  if (i >= n_adapt) break;
37  field ch = navier_stokes_criterion (Re,uh);
38  omega = adapt (ch, options);
39  o.open (omega.name(), "geo");
40  o << omega;
41  o.close();
42  Xh = cavity_space (omega, "P2");
43  Qh = space (omega, "P1");
44  uh = cavity_field (Xh, 1.0);
45  ph = field (Qh, 0);
46  fh = field (Xh, 0);
47  }
48 }
49