/usr/share/doc/rheolef-doc/examples/navier_stokes_taylor_error_dg.cc is in rheolef-doc 6.5-1.
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 | #include "rheolef.h"
using namespace rheolef;
using namespace std;
#include "taylor_exact.icc"
int main(int argc, char**argv) {
environment rheolef(argc, argv);
Float err_u_linf_expected = (argc > 1) ? atof(argv[1]) : 1e+38;
Float err_p_linf_expected = (argc > 2) ? atof(argv[2]) : err_u_linf_expected;
bool have_kinetic_energy = (argc > 3);
bool dump = (argc > 4);
Float Re;
field uh, ph;
din >> catchmark("Re") >> Re
>> catchmark("u") >> uh
>> catchmark("p") >> ph;
space Xh = uh.get_space(),
Qh = ph.get_space();
geo omega = Xh.get_geo();
size_t k = Xh.degree();
size_t d = omega.dimension();
quadrature_option_type qopt;
qopt.set_family(quadrature_option_type::gauss);
qopt.set_order(2*k+1);
#ifdef TODO
Float p_moy = integrate (omega, ph, qopt);
ph = ph-p_moy;
#else // TODO
trial p (Qh); test q (Qh);
form mp = integrate(p*q);
Float p_moy = mp (ph, field(Qh,1));
ph = ph-p_moy;
#endif // TODO
string high_approx = "P"+itos(k+1)+"d";
space Xh1 (omega, high_approx, "vector"),
Qh1 (omega, high_approx);
field euh = interpolate (Xh1, uh-u_exact());
field eph = interpolate (Qh1, ph-p_exact(Re,have_kinetic_energy));
Float err_u_l2 = sqrt(integrate (omega, norm2(uh-u_exact()), qopt));
Float err_u_linf = euh.max_abs();
Float err_u_h1 = sqrt(integrate (omega, norm2(grad_h(euh)), qopt)
+ integrate (omega.sides(), (1/h_local())*norm2(jump(euh)), qopt));
Float err_p_l2 = sqrt(integrate (omega, sqr(ph-p_exact(Re,have_kinetic_energy)), qopt));
Float err_p_linf = eph.max_abs();
derr << "err_u_l2 = " << err_u_l2 << endl
<< "err_u_linf = " << err_u_linf << endl
<< "err_u_h1 = " << err_u_h1 << endl
<< "err_p_l2 = " << err_p_l2 << endl
<< "err_p_linf = " << err_p_linf << endl;
if (dump) {
dout << catchmark("uh") << uh
<< catchmark("u") << interpolate (Xh, u_exact())
<< catchmark("eu") << euh
<< catchmark("ph") << ph
<< catchmark("p") << interpolate (Qh, p_exact(Re,have_kinetic_energy))
<< catchmark("ep") << eph;
}
return ((err_u_linf <= err_u_linf_expected) && (err_p_linf <= err_p_linf_expected)) ? 0 : 1;
}
|