oldroyd_theta_scheme3.h

Go to the documentation of this file.

template <class P>
void oldroyd_theta_scheme<P>::sub_step1 (
  const field& tau_h0, const field& uh0, const field& ph0,
  field& tau_h, field& uh, field& ph) const
{
  update_transport_stress (uh0);
  field gamma_h = inv_mt*(th*tau_h0 - thb);
  test v (Xh), xi (Th);
  field lh = lambda*integrate (dot(uh0,v))
           + b*(c1*tau_h0 + c2*gamma_h);
  ph = ph0;
  uh.set_u() = uh0.u();
  stokes.solve (lh, field(Qh,0), uh, ph);
  field Duh = inv_mt*integrate(ddot(D(uh),xi));
  tau_h = c1*tau_h0 + c2*gamma_h + 2*c3*Duh;
}
template <class P>
void oldroyd_theta_scheme<P>::sub_step2 (
  const field& uh0,
  const field& tau_h1, const field& uh1,
  field& tau_h, field& uh) const
{
  uh = (1-theta)/theta*uh1 - (1-2*theta)/theta*uh0;
  test xi (Th);
  if (We == 0) {
    field Duh = inv_mt*integrate(ddot(D(uh),xi));
    tau_h = 2*alpha*Duh;
    return;
  }
  update_transport_stress (uh);
  form  th_nu = th + nu*mt;
  typename P::tau_upstream tau_up (Th.get_geo(), We, alpha); 
  field lh = integrate (ddot(c4*tau_h1 + 2*c5*D(uh1),xi))
           + integrate ("boundary",
                        max(0, -dot(uh,normal()))*ddot(tau_up,xi));
  problem transport (th_nu);
  transport.solve (lh, tau_h);
}
template <class P>
void
oldroyd_theta_scheme<P>::update_transport_stress (const field& uh) const {
  typename P::tau_upstream tau_up (Th.get_geo(), We, alpha);
  trial tau (Th); test xi (Th);
  auto ma = 0.5*((1-a)*grad(uh) - (1+a)*trans(grad(uh)));
  auto beta_a = tau*ma + trans(ma)*tau;
  th  = integrate (ddot(grad_h(tau)*uh + beta_a,xi))
      + integrate ("boundary", max(0, -dot(uh,normal()))*ddot(tau,xi))
      + integrate ("internal_sides", 
           - dot(uh,normal())*ddot(jump(tau),average(xi))
           + 0.5*abs(dot(uh,normal()))*ddot(jump(tau),jump(xi)));
  thb = integrate ("boundary", max(0, -dot(uh,normal()))*ddot(tau_up,xi));
}