rheolef  6.5
convect.cc
Go to the documentation of this file.
1 #include "rheolef.h"
2 using namespace rheolef;
3 using namespace std;
4 #include "rotating-hill.h"
5 int main (int argc, char **argv) {
6  environment rheolef (argc,argv);
7  geo omega (argv[1]);
8  string approx = (argc > 2) ? argv[2] : "P1";
9  Float nu = (argc > 3) ? atof(argv[3]) : 1e-2;
10  size_t n_max = (argc > 4) ? atoi(argv[4]) : 50;
11  size_t d = omega.dimension();
12  Float delta_t = 2*acos(-1.)/n_max;
13  space Vh (omega, approx, "vector");
14  field uh = interpolate (Vh, u(d));
15  space Xh (omega, approx);
16  Xh.block ("boundary");
17  field phi_h = interpolate (Xh, phi(d,nu,0));
18  characteristic X (-delta_t*uh);
21  qopt.set_order (Xh.degree());
22  trial phi (Xh); test psi (Xh);
23  branch event ("t","phi");
24  dout << catchmark("nu") << nu << endl
25  << event (0, phi_h);
26  for (size_t n = 1; n <= n_max; n++) {
27  Float t = n*delta_t;
28  Float c1 = 1 + delta_t*phi::sigma(d,nu,t);
29  Float c2 = delta_t*nu;
30  form a = integrate (c1*phi*psi + c2*dot(grad(phi),grad(psi)), qopt);
31  field lh = integrate (compose(phi_h, X)*psi, qopt);
32  solver sa (a.uu());
33  phi_h.set_u() = sa.solve (lh.u() - a.ub()*phi_h.b());
34  dout << event (t, phi_h);
35  }
36 }
37