Rheolef  7.1
an efficient C++ finite element environment
burgers_diffusion_dg.cc
Go to the documentation of this file.
1 #include "rheolef.h"
26 using namespace rheolef;
27 using namespace std;
28 #include "burgers.icc"
29 #include "burgers_flux_godunov.icc"
32 #undef NEUMANN
34 int main(int argc, char**argv) {
35  environment rheolef (argc, argv);
36  geo omega (argv[1]);
37  space Xh (omega, argv[2]);
38  size_t k = Xh.degree();
39  Float epsilon = (argc > 3) ? atof(argv[3]) : 0.1;
40  size_t nmax = (argc > 4) ? atoi(argv[4]) : 500;
41  Float tf = (argc > 5) ? atof(argv[5]) : 1;
42  size_t p = (argc > 6) ? atoi(argv[6]) : min(k+1,rk::pmax);
43  Float delta_t = tf/nmax;
44  size_t d = omega.dimension();
45  Float beta = (k+1)*(k+d)/Float(d);
46  trial u (Xh); test v (Xh);
47  form m = integrate (u*v);
48  integrate_option iopt;
49  iopt.invert = true;
50  form inv_m = integrate (u*v, iopt);
51  form a = epsilon*(
52  integrate (dot(grad_h(u),grad_h(v)))
53 #ifdef NEUMANN
54  + integrate ("internal_sides",
55 #else // NEUMANN
56  + integrate ("sides",
57 #endif // NEUMANN
58  beta*penalty()*jump(u)*jump(v)
59  - jump(u)*average(dot(grad_h(v),normal()))
60  - jump(v)*average(dot(grad_h(u),normal()))));
61  vector<problem> pb (p+1);
62  for (size_t i = 1; i <= p; ++i) {
63  form ci = m + delta_t*rk::alpha[p][i][i]*a;
64  pb[i] = problem(ci);
65  }
66  vector<field> uh(p+1, field(Xh,0));
67  uh[0] = interpolate (Xh, u_init(epsilon));
68  branch even("t","u");
69  dout << catchmark("epsilon") << epsilon << endl
70  << even(0,uh[0]);
71  for (size_t n = 0; n < nmax; ++n) {
72  Float tn = n*delta_t;
73  Float t = tn + delta_t;
74  field uh_next = uh[0] - delta_t*rk::tilde_beta[p][0]*(inv_m*gh(epsilon, tn, uh[0], v));
75  for (size_t i = 1; i <= p; ++i) {
76  Float ti = tn + rk::gamma[p][i]*delta_t;
77  field rhs = m*uh[0] - delta_t*rk::tilde_alpha[p][i][0]*gh(epsilon, tn, uh[0], v);
78  for (size_t j = 1; j <= i-1; ++j) {
79  Float tj = tn + rk::gamma[p][j]*delta_t;
80  rhs -= delta_t*( rk::alpha[p][i][j]*(a*uh[j] - lh(epsilon,tj,v))
81  + rk::tilde_alpha[p][i][j]*gh(epsilon, tj, uh[j], v));
82  }
83  rhs += delta_t*rk::alpha[p][i][i]*lh (epsilon, ti, v);
84  pb[i].solve (rhs, uh[i]);
85  uh_next -= delta_t*(inv_m*( rk::beta[p][i]*(a*uh[i] - lh(epsilon,ti,v))
86  + rk::tilde_beta[p][i]*gh(epsilon, ti, uh[i], v)));
87  }
88  uh_next = limiter(uh_next);
89  dout << even(tn+delta_t,uh_next);
90  uh[0] = uh_next;
91  }
92 }
The Burgers equation – the f function.
int main(int argc, char **argv)
The diffusive Burgers equation – its exact solution.
u_exact u_init
The diffusive Burgers equation – operators.
field lh(Float epsilon, Float t, const test &v)
field gh(Float epsilon, Float t, const field &uh, const test &v)
The Burgers equation – the Godonov flux.
see the Float page for the full documentation
see the branch page for the full documentation
see the field page for the full documentation
see the form page for the full documentation
see the geo page for the full documentation
see the catchmark page for the full documentation
Definition: catchmark.h:67
see the environment page for the full documentation
Definition: environment.h:104
see the integrate_option page for the full documentation
problem_basic< Float > problem
Definition: problem.h:163
double Float
see the Float page for the full documentation
Definition: Float.h:143
odiststream dout(cout)
see the diststream page for the full documentation
Definition: diststream.h:430
field_basic< Float > field
see the field page for the full documentation
Definition: field.h:419
see the space page for the full documentation
see the test page for the full documentation
see the test page for the full documentation
This file is part of Rheolef.
field_basic< T, M > limiter(const field_basic< T, M > &uh, const T &bar_g_S, const limiter_option &opt)
see the limiter page for the full documentation
Definition: limiter.cc:65
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
Definition: vec_expr_v2.h:415
std::enable_if< details::is_field_convertible< Expr >::value,details::field_expr_v2_nonlinear_terminal_field< typename Expr::scalar_type,typename Expr::memory_type,details::differentiate_option::gradient >>::type grad_h(const Expr &expr)
grad_h(uh): see the expression page for the full documentation
std::enable_if< details::is_field_expr_v2_nonlinear_arg< Expr >::value &&! is_undeterminated< Result >::value, Result >::type integrate(const geo_basic< T, M > &omega, const Expr &expr, const integrate_option &iopt, Result dummy=Result())
see the integrate page for the full documentation
Definition: integrate.h:202
field_basic< T, M > interpolate(const space_basic< T, M > &V2h, const field_basic< T, M > &u1h)
see the interpolate page for the full documentation
Definition: interpolate.cc:233
details::field_expr_v2_nonlinear_terminal_function< details::penalty_pseudo_function< Float > > penalty()
penalty(): see the expression page for the full documentation
details::field_expr_v2_nonlinear_terminal_function< details::normal_pseudo_function< Float > > normal()
normal: see the expression page for the full documentation
Float tilde_alpha[][pmax+1][pmax+1]
Float tilde_beta[][pmax+1]
Float gamma[][pmax+1]
Float beta[][pmax+1]
Float alpha[][pmax+1][pmax+1]
constexpr size_t pmax
rheolef - reference manual
The semi-implicit Runge-Kutta scheme – coefficients.
Definition: sphere.icc:25
Definition: leveque.h:25
Float u(const point &x)
Float epsilon