burgers_diffusion_dg.cc

The diffusive Burgers equation by the discontinuous Galerkin method

#include "rheolef.h"
using namespace rheolef;
using namespace std;
#include "burgers.icc"
#include "burgers_flux_godunov.icc"
#include "runge_kutta_semiimplicit.icc"
#include "burgers_diffusion_exact.h"
#undef NEUMANN
#include "burgers_diffusion_operators.icc"
int main(int argc, char**argv) {
environment rheolef (argc, argv);
geo omega (argv[1]);
space Xh (omega, argv[2]);
size_t k = Xh.degree();
Float epsilon = (argc > 3) ? atof(argv[3]) : 0.1;
size_t nmax = (argc > 4) ? atoi(argv[4]) : 500;
Float tf = (argc > 5) ? atof(argv[5]) : 1;
size_t p = (argc > 6) ? atoi(argv[6]) : min(k+1,rk::pmax);
Float delta_t = tf/nmax;
size_t d = omega.dimension();
Float beta = (k+1)*(k+d)/Float(d);
trial u (Xh); test v (Xh);
form m = integrate (u*v);
integrate_option iopt;
iopt.invert = true;
form inv_m = integrate (u*v, iopt);
form a = epsilon*(
integrate (dot(grad_h(u),grad_h(v)))
#ifdef NEUMANN
+ integrate ("internal_sides",
#else // NEUMANN
+ integrate ("sides",
#endif // NEUMANN
beta*penalty()*jump(u)*jump(v)
- jump(u)*average(dot(grad_h(v),normal()))
- jump(v)*average(dot(grad_h(u),normal()))));
vector<problem> pb (p+1);
for (size_t i = 1; i <= p; ++i) {
form ci = m + delta_t*rk::alpha[p][i][i]*a;
pb[i] = problem(ci);
}
vector<field> uh(p+1, field(Xh,0));
uh[0] = interpolate (Xh, u_init(epsilon));
branch even("t","u");
dout << catchmark("epsilon") << epsilon << endl
<< even(0,uh[0]);
for (size_t n = 0; n < nmax; ++n) {
Float tn = n*delta_t;
Float t = tn + delta_t;
field uh_next = uh[0] - delta_t*rk::tilde_beta[p][0]*(inv_m*gh(epsilon, tn, uh[0], v));
for (size_t i = 1; i <= p; ++i) {
Float ti = tn + rk::gamma[p][i]*delta_t;
field rhs = m*uh[0] - delta_t*rk::tilde_alpha[p][i][0]*gh(epsilon, tn, uh[0], v);
for (size_t j = 1; j <= i-1; ++j) {
Float tj = tn + rk::gamma[p][j]*delta_t;
rhs -= delta_t*( rk::alpha[p][i][j]*(a*uh[j] - lh(epsilon,tj,v))
+ rk::tilde_alpha[p][i][j]*gh(epsilon, tj, uh[j], v));
}
rhs += delta_t*rk::alpha[p][i][i]*lh (epsilon, ti, v);
pb[i].solve (rhs, uh[i]);
uh_next -= delta_t*(inv_m*( rk::beta[p][i]*(a*uh[i] - lh(epsilon,ti,v))
+ rk::tilde_beta[p][i]*gh(epsilon, ti, uh[i], v)));
}
uh_next = limiter(uh_next);
dout << even(tn+delta_t,uh_next);
uh[0] = uh_next;
}
}
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