Loading...
Searching...
No Matches
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();
}
gh
field gh(Float epsilon, Float t, const field &uh, const test &v)
Definition burgers_diffusion_operators.icc:37
Float
see the Float page for the full documentation
field
see the field page for the full documentation
form
see the form page for the full documentation
geo
see the geo page for the full documentation
problem_mixed
see the problem_mixed page for the full documentation
test
see the test page for the full documentation
trial
see the test page for the full documentation
u
point u(const point &x)
Definition diffusion_transport_tensor_exact.icc:25
residue
field residue(Float p, const field &uh)
Definition p_laplacian_post.cc:35
oldroyd_theta_scheme::inv_mt
form inv_mt
Definition oldroyd_theta_scheme.h:46
oldroyd_theta_scheme::th
form th
Definition oldroyd_theta_scheme.h:47
oldroyd_theta_scheme::mp
form mp
Definition oldroyd_theta_scheme.h:46
oldroyd_theta_scheme::We
Float We
Definition oldroyd_theta_scheme.h:42
oldroyd_theta_scheme::delta_t
Float delta_t
Definition oldroyd_theta_scheme.h:42
oldroyd_theta_scheme::Xh
space Xh
Definition oldroyd_theta_scheme.h:45
oldroyd_theta_scheme::sub_step1
void sub_step1(const field &tau_h0, const field &uh0, const field &ph0, field &tau_h, field &uh, field &ph) const
Definition oldroyd_theta_scheme3.h:26
oldroyd_theta_scheme::thb
field thb
Definition oldroyd_theta_scheme.h:48
oldroyd_theta_scheme::alpha
Float alpha
Definition oldroyd_theta_scheme.h:42
oldroyd_theta_scheme::d
form d
Definition oldroyd_theta_scheme.h:46
oldroyd_theta_scheme::tol
Float tol
Definition oldroyd_theta_scheme.h:42
oldroyd_theta_scheme::b
form b
Definition oldroyd_theta_scheme.h:46
oldroyd_theta_scheme::solve
bool solve(field &tau_h, field &uh, field &ph)
Definition oldroyd_theta_scheme2.h:26
oldroyd_theta_scheme::initial
void initial(const geo &omega, field &tau_h, field &uh, field &ph, string restart)
Definition oldroyd_theta_scheme2.h:47
oldroyd_theta_scheme::Th
space Th
Definition oldroyd_theta_scheme.h:45
oldroyd_theta_scheme::residue
Float residue(field &tau_h, field &uh, field &ph) const
Definition oldroyd_theta_scheme2.h:81
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::max_iter
size_t max_iter
Definition oldroyd_theta_scheme.h:43
oldroyd_theta_scheme::reset
void reset(const geo &omega)
Definition oldroyd_theta_scheme1.h:31
oldroyd_theta_scheme::Qh
space Qh
Definition oldroyd_theta_scheme.h:45
oldroyd_theta_scheme::update_transport_stress
void update_transport_stress(const field &uh) const
Definition oldroyd_theta_scheme3.h:65
oldroyd_theta_scheme::sub_step2
void sub_step2(const field &uh0, const field &tau_h1, const field &uh1, field &tau_h, field &uh) const
Definition oldroyd_theta_scheme3.h:42
u
Definition leveque.h:25