The Oldroyd problem by the theta-scheme – class body.
template <class P>
reset (uh.get_geo());
field tau_h0 = tau_h, uh0 = uh, ph0 = ph;
derr << "# n t rel_err residue lambda_min" << endl;
derr << "0 0 0 " << r << endl;
for (size_t n = 1; n <= max_iter; ++n) {
step (tau_h0, uh0, ph0, tau_h, uh, ph);
Float rel_err_prec = rel_err, r_prec = r;
rel_err =
field(tau_h-tau_h0).max_abs() +
field(uh-uh0).max_abs();
derr << n << " " << n*delta_t << " " << rel_err << " " << r << endl;
if (rel_err < tol) return true;
if (rel_err_prec != 0 && ((rel_err > 10*rel_err_prec && r > 10*r_prec) ||
(rel_err > 1e5 && r > 1e5) )) return false;
tau_h0 = tau_h; uh0 = uh; ph0 = ph;
}
return (rel_err < sqrt(tol));
}
template <class P>
reset (omega);
if (restart == "") {
uh = P::velocity_field (Xh);
form c0 = integrate (2*ddot(D(
u),D(v)));
s0.set_metric(mp);
field Duh = inv_mt*integrate(ddot(D(uh),xi));
tau_h = 2*alpha*Duh;
} else {
idiststream in (restart, "field");
in >> catchmark("tau") >> tau_h
>> catchmark("u") >> uh
>> catchmark("p") >> ph;
}
}
template <class P>
field tau_h1 = tau_h0, uh1 = uh0, ph1 = ph0;
sub_step1 (tau_h0, uh0, ph0, tau_h1, uh1, ph1);
field tau_h2 = tau_h1, uh2 = uh1;
sub_step2 (uh0, tau_h1, uh1, tau_h2, uh2);
sub_step1 (tau_h2, uh2, ph1, tau_h, uh, ph);
}
template <class P>
update_transport_stress (uh);
field rt = We*(th*tau_h-thb) + integrate (
ddot(tau_h - alpha*
gh, xi));
ru.set_b() = 0;
return rt.u().max_abs() + ru.u().max_abs() + rp.u().max_abs();
}
field gh(Float epsilon, Float t, const field &uh, const test &v)
see the Float page for the full documentation
see the field page for the full documentation
see the geo page for the full documentation
see the problem_mixed page for the full documentation
see the test page for the full documentation
see the test page for the full documentation
Float alpha[pmax+1][pmax+1]
std::enable_if< details::has_field_rdof_interface< Expr >::value, details::field_expr_v2_nonlinear_terminal_field< typenameExpr::scalar_type, typenameExpr::memory_type, details::differentiate_option::gradient > >::type D(const Expr &expr)
D(uh): see the expression page for the full documentation.
T ddot(const tensor_basic< T > &a, const tensor_basic< T > &b)
ddot(x,y): see the expression page for the full documentation
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
field residue(Float p, const field &uh)
bool solve(field &tau_h, field &uh, field &ph)
void initial(const geo &omega, field &tau_h, field &uh, field &ph, string restart)
Float residue(field &tau_h, field &uh, field &ph) const
void step(const field &tau_h0, const field &uh0, const field &ph0, field &tau_h, field &uh, field &ph) const