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 }
rheolef::problem
problem_basic< Float > problem
Definition: problem.h:163
burgers_flux_godunov.icc
The Burgers equation – the Godonov flux.
rk::pmax
constexpr size_t pmax
Definition: runge_kutta_semiimplicit.icc:27
form
see the form page for the full documentation
gh
field gh(Float epsilon, Float t, const field &uh, const test &v)
Definition: burgers_diffusion_operators.icc:37
rheolef::catchmark
see the catchmark page for the full documentation
Definition: catchmark.h:67
rheolef::dot
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
field
see the field page for the full documentation
rheolef::normal
details::field_expr_v2_nonlinear_terminal_function< details::normal_pseudo_function< Float > > normal()
normal: see the expression page for the full documentation
Definition: field_expr_terminal.h:439
rheolef::field
field_basic< Float > field
see the field page for the full documentation
Definition: field.h:419
rheolef::integrate
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
u_init
u_exact u_init
Definition: burgers_diffusion_exact.h:32
space
see the space page for the full documentation
rk::tilde_alpha
Float tilde_alpha[][pmax+1][pmax+1]
Definition: runge_kutta_semiimplicit.icc:48
rk::gamma
Float gamma[][pmax+1]
Definition: runge_kutta_semiimplicit.icc:70
rheolef.h
rheolef - reference manual
p
Definition: sphere.icc:25
rheolef::integrate_option
see the integrate_option page for the full documentation
Definition: integrate_option.h:125
burgers_diffusion_operators.icc
The diffusive Burgers equation – operators.
a
Definition: diffusion_isotropic.h:25
rheolef::interpolate
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
rheolef::environment
see the environment page for the full documentation
Definition: environment.h:115
lh
field lh(Float epsilon, Float t, const test &v)
Definition: burgers_diffusion_operators.icc:25
mkgeo_sector.m
m
Definition: mkgeo_sector.sh:118
rheolef::limiter
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
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::integrate_option::invert
bool invert
Definition: integrate_option.h:168
test
see the test page for the full documentation
u
Definition: leveque.h:25
Float
see the Float page for the full documentation
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
branch
see the branch page for the full documentation
u
Float u(const point &x)
Definition: transmission_error.cc:26
burgers.icc
The Burgers equation – the f function.
mkgeo_ball.a
a
Definition: mkgeo_ball.sh:151
rheolef::penalty
details::field_expr_v2_nonlinear_terminal_function< details::penalty_pseudo_function< Float > > penalty()
penalty(): see the expression page for the full documentation
Definition: field_expr_terminal.h:626
burgers_diffusion_exact.h
The diffusive Burgers equation – its exact solution.
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
runge_kutta_semiimplicit.icc
The semi-implicit Runge-Kutta scheme – coefficients.
rk::tilde_beta
Float tilde_beta[][pmax+1]
Definition: runge_kutta_semiimplicit.icc:65
rheolef::grad_h
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
Definition: field_expr_terminal.h:949
trial
see the test page for the full documentation
rheolef::Float
double Float
see the Float page for the full documentation
Definition: Float.h:143
rk::alpha
Float alpha[][pmax+1][pmax+1]
Definition: runge_kutta_semiimplicit.icc:36
rheolef::dout
odiststream dout(cout)
see the diststream page for the full documentation
Definition: diststream.h:430
epsilon
Float epsilon
Definition: transmission_error.cc:25
rheolef::std
Definition: vec_expr_v2.h:391
geo
see the geo page for the full documentation
rk::beta
Float beta[][pmax+1]
Definition: runge_kutta_semiimplicit.icc:60
main
int main(int argc, char **argv)
Definition: burgers_diffusion_dg.cc:34