/usr/share/doc/rheolef-doc/examples/burgers_diffusion_dg.cc is in rheolef-doc 6.7-6.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 | #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.icc"
#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)/d;
trial u (Xh); test v (Xh);
form m = integrate (u*v);
form_option_type fopt;
fopt.invert = true;
form inv_m = integrate (u*v, fopt);
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<solver> sc (p+1);
for (size_t i = 1; i <= p; ++i) {
form ci = m + delta_t*rk::alpha[p][i][i]*a;
sc[i] = solver(ci.uu());
}
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);
uh[i].set_u() = sc[i].solve (rhs.u());
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;
}
}
|