Rheolef  7.1
an efficient C++ finite element environment
oldroyd_theta_scheme2.h
Go to the documentation of this file.
1 template <class P>
27  reset (uh.get_geo());
28  field tau_h0 = tau_h, uh0 = uh, ph0 = ph;
29  derr << "# n t rel_err residue lambda_min" << endl;
30  Float r = residue (tau_h, uh, ph);
31  Float rel_err = 0;
32  derr << "0 0 0 " << r << endl;
33  for (size_t n = 1; n <= max_iter; ++n) {
34  step (tau_h0, uh0, ph0, tau_h, uh, ph);
35  Float rel_err_prec = rel_err, r_prec = r;
36  r = residue (tau_h, uh, ph);
37  rel_err = field(tau_h-tau_h0).max_abs() + field(uh-uh0).max_abs();
38  derr << n << " " << n*delta_t << " " << rel_err << " " << r << endl;
39  if (rel_err < tol) return true;
40  if (rel_err_prec != 0 && ((rel_err > 10*rel_err_prec && r > 10*r_prec) ||
41  (rel_err > 1e5 && r > 1e5) )) return false;
42  tau_h0 = tau_h; uh0 = uh; ph0 = ph;
43  }
44  return (rel_err < sqrt(tol));
45 }
46 template <class P>
48  const geo& omega, field& tau_h, field& uh, field& ph, string restart) {
49  reset (omega);
50  ph = field(Qh,0);
51  if (restart == "") {
52  uh = P::velocity_field (Xh);
53  trial u (Xh); test v (Xh), xi(Th);
54  form c0 = integrate (2*ddot(D(u),D(v)));
55  problem_mixed s0 (c0, d);
56  s0.set_metric(mp);
57  s0.solve (field(Xh,0), field(Qh,0), uh, ph);
58  field Duh = inv_mt*integrate(ddot(D(uh),xi));
59  tau_h = 2*alpha*Duh;
60  } else {
61  tau_h = field(Th);
62  uh = field(Xh);
63  idiststream in (restart, "field");
64  in >> catchmark("tau") >> tau_h
65  >> catchmark("u") >> uh
66  >> catchmark("p") >> ph;
67  }
68 }
69 template <class P>
71  const field& tau_h0, const field& uh0, const field& ph0,
72  field& tau_h, field& uh, field& ph) const {
73  field tau_h1 = tau_h0, uh1 = uh0, ph1 = ph0;
74  sub_step1 (tau_h0, uh0, ph0, tau_h1, uh1, ph1);
75  field tau_h2 = tau_h1, uh2 = uh1;
76  sub_step2 (uh0, tau_h1, uh1, tau_h2, uh2);
77  sub_step1 (tau_h2, uh2, ph1, tau_h, uh, ph);
78 }
79 template <class P>
80 Float
82  update_transport_stress (uh);
83  test xi (Th);
84  field Duh = inv_mt*integrate(ddot(D(uh),xi));
85  field gh = 2*Duh;
86  field rt = We*(th*tau_h-thb) + integrate (ddot(tau_h - alpha*gh, xi));
87  field ru = b*(tau_h + (1-alpha)*gh) - d.trans_mult(ph);
88  ru.set_b() = 0;
89  field rp = d*uh;
90  return rt.u().max_abs() + ru.u().max_abs() + rp.u().max_abs();
91 }
mkgeo_ball.n
int n
Definition: mkgeo_ball.sh:150
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
bdf::alpha
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
field
see the field page for the full documentation
rheolef::details::reset
void reset(T &x)
Definition: keller_details.h:31
residue
field residue(Float p, const field &uh)
Definition: p_laplacian_post.cc:35
oldroyd_theta_scheme::solve
bool solve(field &tau_h, field &uh, field &ph)
Definition: oldroyd_theta_scheme2.h:26
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
problem_mixed
see the problem_mixed page for the full documentation
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
oldroyd_theta_scheme::step
void step(const field &tau_h0, const field &uh0, const field &ph0, field &tau_h, field &uh, field &ph) const
Definition: oldroyd_theta_scheme2.h:70
oldroyd_theta_scheme::initial
void initial(const geo &omega, field &tau_h, field &uh, field &ph, string restart)
Definition: oldroyd_theta_scheme2.h:47
rheolef::ddot
T ddot(const tensor_basic< T > &a, const tensor_basic< T > &b)
ddot(x,y): see the expression page for the full documentation
Definition: tensor.cc:278
test
see the test page for the full documentation
u
Definition: leveque.h:25
Float
see the Float page for the full documentation
rheolef::D
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 D(const Expr &expr)
D(uh): see the expression page for the full documentation.
Definition: field_expr_terminal.h:969
u
Float u(const point &x)
Definition: transmission_error.cc:26
trial
see the test page for the full documentation
oldroyd_theta_scheme::residue
Float residue(field &tau_h, field &uh, field &ph) const
Definition: oldroyd_theta_scheme2.h:81
geo
see the geo page for the full documentation