rheolef  6.5
convect2.cc
Go to the documentation of this file.
1 // mkgeo_grid -e 1000 -a -2 -b 2 > line2.geo
2 // ./convect line2 1e-3 5 > line2.branch
3 // ./convect2 line2 1e-3 5 > line22.branch
4 // (./convect_error < line2.branch > line2-cmp.branch) 2>&1|g -v trace |g -v load |tee cv.gdat
5 // (./convect_error < line22.branch > line22-cmp.branch) 2>&1|g -v trace |g -v load |tee cv2.gdat
6 #include "rheolef.h"
7 using namespace rheolef;
8 using namespace std;
9 #include "rotating-hill.h"
10 int main (int argc, char **argv) {
11  environment rheolef (argc,argv);
12  geo omega (argv[1]);
13  string approx = (argc > 2) ? argv[2] : "P1";
14  Float nu = (argc > 3) ? atof(argv[3]) : 1e-2;
15  size_t n_max = (argc > 4) ? atoi(argv[4]) : 50;
16  size_t d = omega.dimension();
17  Float delta_t = 2*acos(-1.)/n_max;
18  space Vh (omega, approx, "vector");
19  field uh = interpolate (Vh, u(d));
20  space Xh (omega, approx);
21  Xh.block ("boundary");
22  field phi_h = interpolate (Xh, phi(d,nu,0));
23  field phi_h_prec = interpolate (Xh, phi(d,nu,-delta_t));
24  characteristic X1 ( -delta_t*uh);
25  characteristic X2 (-2*delta_t*uh);
28  qopt.set_order (Xh.degree());
29  trial varphi (Xh); test psi (Xh);
30  branch event ("t","phi");
31  dout << catchmark("nu") << nu << endl
32  << event (0, phi_h);
33  for (size_t n = 1; n <= n_max; n++) {
34  Float t = n*delta_t;
35  Float alpha = 1.5 + delta_t*phi::sigma(d,nu,t);
36  Float beta = delta_t*nu;
37  form c = integrate (alpha*varphi*psi + beta*dot(grad(varphi),grad(psi)), qopt);
38  field lh = 2.0*integrate (compose(phi_h, X1)*psi, qopt)
39  - 0.5*integrate (compose(phi_h_prec, X2)*psi, qopt);
40 #ifdef TODO // boost::proto & g++-4.8
41  field lh = integrate ((2.*compose(phi_h, X1) - 0.5*compose(phi_h_prec, X2))*psi, qopt);
42 #endif // TODO
43  solver sc (c.uu());
44  phi_h_prec = phi_h;
45  phi_h.set_u() = sc.solve (lh.u() - c.ub()*phi_h.b());
46  dout << event (t, phi_h);
47  }
48 }
49