oldroyd_theta_scheme2.h

Go to the documentation of this file.

template <class P>
bool oldroyd_theta_scheme<P>::solve(field& tau_h, field& uh, field& ph) {
  reset (uh.get_geo());
  field tau_h0 = tau_h, uh0 = uh, ph0 = ph;
  derr << "# n t rel_err residue lambda_min" << endl;
  Float r = residue (tau_h, uh, ph);
  Float rel_err = 0;
  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;
    r = residue (tau_h, uh, ph);
    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>
void oldroyd_theta_scheme<P>::initial (
    const geo& omega, field& tau_h, field& uh, field& ph, string restart) {
  reset (omega);
  ph = field(Qh,0);
  if (restart == "") {
    uh = P::velocity_field (Xh);
    trial u (Xh); test v (Xh), xi(Th);
    form c0  = integrate (2*ddot(D(u),D(v)));
    problem_mixed s0 (c0, d);
    s0.set_metric(mp);
    s0.solve (field(Xh,0), field(Qh,0), uh, ph); 
    field Duh = inv_mt*integrate(ddot(D(uh),xi));
    tau_h = 2*alpha*Duh;
  } else {
    tau_h = field(Th);
    uh    = field(Xh);
    idiststream in (restart, "field");
    in >> catchmark("tau") >> tau_h
       >> catchmark("u")   >> uh
       >> catchmark("p")   >> ph;
  }
}
template <class P>
void oldroyd_theta_scheme<P>::step (
    const field& tau_h0, const field& uh0, const field& ph0,
    field& tau_h, field& uh, field& ph) const {
  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>
Float
oldroyd_theta_scheme<P>::residue(field& tau_h, field& uh, field& ph) const{
  update_transport_stress (uh);
  test xi (Th);
  field Duh = inv_mt*integrate(ddot(D(uh),xi));
  field gh = 2*Duh;
  field rt = We*(th*tau_h-thb) + integrate (ddot(tau_h - alpha*gh, xi));
  field ru = b*(tau_h + (1-alpha)*gh) - d.trans_mult(ph);
  ru.set_b() = 0;
  field rp = d*uh;
  return rt.u().max_abs() + ru.u().max_abs() + rp.u().max_abs();
}