rheolef  6.6
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,0));
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 = Float(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 = Float(2.0)*integrate (compose(phi_h, X1)*psi, qopt)
39  - Float(0.5)*integrate (compose(phi_h_prec, X2)*psi, qopt);
40 #ifdef TODO // boost::proto & g++-4.8
41  field lh = integrate ((Float(2)*compose(phi_h, X1) - Float(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 }
quadrature_option_type - send options to the integrate function
boost::proto::result_of::make_expr< boost::proto::tag::function, vec_domain, const vec_detail::compose_< Function >, const vec< T, M > & >::type compose(const Function &f, const vec< T, M > &x)
field - piecewise polynomial finite element field
Definition: field.h:225
const vec< T, M > & u() const
Definition: field.h:279
field_vf_expr< field_vf_expr_grad< test_basic< T, M, VfTag > >> grad(const test_basic< T, M, VfTag > &x)
Definition: operators2.h:74
STL namespace.
vec< T, M > & set_u()
Definition: field.h:281
irheostream, orheostream - large data streams
Definition: compiler.h:7
const csr< T, M > & uu() const
Definition: form.h:140
static field_basic< T, sequential > interpolate(const space_basic< T, sequential > &Vh, const field_basic< T, sequential > &uh)
branch - a parameter-dependent sequence of field
Definition: branch.h:46
const csr< T, M > & ub() const
Definition: form.h:141
space – piecewise polynomial finite element space
Definition: space.h:229
double Float
Definition: compiler.h:177
size_t d
field_nonlinear_expr< field_nonlinear_expr_uf< details::acos_,field_expr_terminal_field< T, M > > > acos(const field_basic< T, M > &x)
const vec< T, M > & b() const
Definition: field.h:280
catchmark - iostream manipulator
Definition: catchmark.h:30
Float sigma
T dot(const vec< T, M > &x, const int &y)
form - representation of a finite element bilinear form
Definition: form.h:98
solver - direct or interative solver interface
Definition: solver.h:8
characteristic - the Lagrange-Galerkin method implemented
int main(int argc, char **argv)
Definition: convect2.cc:10
Float u(const point &x)
odiststream dout(cout)
Definition: diststream.h:317
field_basic< T, M > integrate(const geo_basic< T, M > &domain, const test_basic< T, M, details::vf_tag_01 > &expr, const quadrature_option_type &qopt=quadrature_option_type())
integrate - integrate a function or an expression
Definition: integrate.h:90