/usr/share/doc/rheolef-doc/examples/navier_stokes_solve.icc 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 | using namespace std;
int navier_stokes_solve (
Float Re, Float delta_t, field l0h, field& uh, field& ph,
size_t& max_iter, Float& tol, odiststream *p_derr=0) {
const space& Xh = uh.get_space();
const space& Qh = ph.get_space();
string label = "navier-stokes-" + Xh.get_geo().name();
quadrature_option_type qopt;
qopt.set_family(quadrature_option_type::gauss_lobatto);
qopt.set_order(Xh.degree());
trial u (Xh), p (Qh);
test v (Xh), q (Qh);
form mp = integrate (p*q, qopt);
form m = integrate (dot(u,v), qopt);
form a = integrate (2*ddot(D(u),D(v)) + 1.5*(Re/delta_t)*dot(u,v), qopt);
form b = integrate (-div(u)*q, qopt);
solver sa (a.uu());
solver_abtb stokes (a.uu(), b.uu(), mp.uu());
if (p_derr != 0) *p_derr << "[" << label << "] #n |du/dt|" << endl;
field uh1 = uh;
for (size_t n = 0; true; n++) {
field uh2 = uh1;
uh1 = uh;
field uh_star = 2.0*uh1 - uh2;
characteristic X1 ( -delta_t*uh_star);
characteristic X2 (-2.0*delta_t*uh_star);
field l1h = integrate (dot(compose(uh1,X1),v), qopt);
field l2h = integrate (dot(compose(uh2,X2),v), qopt);
field lh = l0h + (Re/delta_t)*(2*l1h - 0.5*l2h);
stokes.solve (lh.u() - a.ub()*uh.b(), -(b.ub()*uh.b()),
uh.set_u(), ph.set_u());
field duh_dt = (3*uh - 4*uh1 + uh2)/(2*delta_t);
Float residual = sqrt(m(duh_dt,duh_dt));
if (p_derr != 0) *p_derr << "[" << label << "] "<< n << " " << residual << endl;
if (residual < tol) {
tol = residual;
max_iter = n;
return 0;
}
if (n == max_iter-1) {
tol = residual;
return 1;
}
}
}
|